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ABSTRACT 


The  numerical  method  developed  by  Chen  and  Patel  (IIHR  Report  No.  285, 
April  1985)  for  the  solution  of  the  partial ly-parabolic  Reynolds-averaged 
Navier-Stokes  equations  has  been  generalized  to  solve  the  fully-elliptic 
equations.  This  method  is  applied  to  calculate  the  flow  over  the  stern  and  in 
the  wake  of  several  ship  forms  for  which  extensive  data  are  available.  The 
report  also  provides  an  overview  of  the  present  status  of  experiments  and 
computational  capabilities  for  such  flows. 
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NOMENCLATURE 


A  grid-attraction  amplitude  function  defined  in  equation  (41) 

A,B,C,R  convection  coefficients  in  the  linearized  convective-transport 

equation  (45) 

a,b,c  (1)  modified  grid-control  functions  defined  by  equation  (36) 

(2)  constants  in  analytical  solution  of  the  one-dimensional 
equations  (49)  and  (51) 

convection  coefficients  in  the  linearized  convective-transport 
equation  for  $(=U,V,W,k,e) 

ap,  ad>  etc.  finite-analytic  coefficients  for  pressure  and  pressure-correc¬ 
tion  equations 

b^.c^.d^  constants  in  transport  equations  (16)  for  4>(=U,V,W,k,e ) 

b0  two-dimensional  grid-control  function  defined  by  equation  (40) 

bc  three-dimensional  grid  correction  function  defined  by  equation 

(41) 

b^  geometric  coefficients  as  defined  in  equation  (27) 

C  dimensionless  crossflow  velocity,  normalized  by  U0 

c  grid-adjustment  factor  defined  in  equation  (4jl.) 

Cb  block  coefficient 

Cnb  finite  analytic  coefficients  for  transport  ‘equations 

(nb  =  NE,  NW,  SE,  SW,  EC,  WC,  NC,-  SC) 

Cd»Cn,Ce  finite-analytic  coefficients  for  pressure  and  pressure-correc¬ 

tion  equations 

Cf  pUo»  coefficient  (=2U^) 

1  2 

Cp  =2p,  pressure  coefficient  normalized  by  j  pUq 


iv 


Cp»C(j,CD 


Cu’Ce1’Ce2 

d 

D1 

A 

D 

* 

D 


jd,dn’ 


etc 


f1 

G 

g 


L 

1 

N 

N 


finite  analytic  coefficients  for  transport  equations  at  cen 
tral,  upstream  and  downstream  nodes 
turbulence-model  constants 
constant  defined  in  (41) 

partial  mass  source  term  defined  in  equation  (60a) 
source  function  in  pressure  equation  (61) 
mass  source  term  in  pressure-correction  equation  (65) 
finite-analytic  coefficients  for  pressure  and  pressure-correc¬ 
tion  equations 

series  summation  term  in  equation  (52) 
grid-control  functions  in  equations  (19)  and  (20) 
turbulence  generation  term  in  equation  (9)  and  (26f) 

(1)  geometric  coefficient  in  equation  (29) 

(2)  source  function  in  linearized  transport  equations  (50)  and 
(52) 

metric  tensor  in  general  curvilinear  coordinates 
conjugate  metric  tensor  in  general  curvilinear  coordinates 
grid  size  in  equation  (46) 

metric  coefficients  or  scale  factors  in  orthogonal  coordinates 
x1 

Jacobian 

2 

(1)  dimensionless  turbulent  kinetic- energy,  normalized  by  UQ 

(2)  grid  size  in  equation  (46) 
length  scale  (body  length) 
grid  size  in  equation  (46) 

distance  normal  to  the  body  cross-section 
UqN/v 


v 


j]1  unit  vectors  normal  to  the  coordinate  lines 

2 

p  dimensionless  pressure,  normalized  by  pUq 

★ 

p  guessed  (imperfect)  pressure  field 

* 

p  p-p  ,  pressure-correction 

3  dimensionless  total  velocity  vector,  normalized  by  U0 

* 

q  projection  of  the  velocity  vector  on  the  body  surface 

/•v 

q  magnitude  of  q 

/V 

q  ,qn,qC  dimensionless  component  velocities  along  £,n,c  directions, 

respectively 
R  radius 

r  ff/L,  dimensionless  radius 


r 

rs»rmax 

Re 


t 


U,V,W 

'k  ic  ★ 

u  ,v  ,w 
u,v,w 

AAA 

u,v,w 

uo 

uc 

u 

IV 

Ux 

uu,  uv,  etc. 


position  vector 

dimensionless  radius  of  the  body  surface  and  outer  boundaries 

UqL/v,  Reynolds  nimber 

1/Re  +  v Jo  as  in  equation  (10) 

t  9 

source  functions  for  transport  quantities  $(=  U,V,W,k,e) 

dimensionless  time,  normalized  by  L/UQ 

dimensionless  velocity  components,  normalized  by  U0 

* 

velocities  obtained  from  guessed  pressure  field  p 
pseudovelocities  in  equation  (57) 
modified  pseudovelocities  in  equation  (59) 
constant  free-streram  (reference)  velocity 
wake  centerline  velocity  normalized  by  U0 
velocity  vector 

(T  /pin)  ,  normalized  friction  (wall  shear)  velocity 

W  U 

? 

dimensionless  Reynolds  stresses,  normalized  by  UQ 
velocity  at  the  edge  of  boundary  layer 


vi 


*,Y,7 

x,y,z;X,Y,Z 

x1 

x,r,8 

y+ 

* 

y 


Cartesian  coordinates 

dimensionless  Cartesian  coordinates 

dimension'ess  orthogonal  coordinates  ( i=l ,2,3) 

dimensionless  cyl indrical-polar  coordinates 

ReUTy,  dimensionless  distance  normal  to  the  wall 

dimensionless  distance  from  the  hull  surface  in  the  y-direction 


Greek 


a 

i 


6 

Y 

5 

e 

K 

A 

m 

v 


*  *  * 
5  ,n  ,5 


7T 


4> 


angle  between  and  nn 

coefficients  in  transport  equation  (22)  for  4>(=U,V,W,k,e )  and 
i=l  ,3,3 

* 

angle  between  q  and  5-direction  on  the  body  surface 
angle  between  t  and  t 
boundary-layer  thickness 

3 

rate  of  turbulent  energy  dissipation,  normalized  by  UQ/L 
von  Karman  constant 
eigenvalues 
kinematic  viscosity 

turbulent  eddy-viscosity,  normalized  by  UQL 

transformed  (general  curvilinear)  coordinates 

general  curvilinear  coordinates  ( i=l ,2,3) 

transformed  coordinates  in  linearized  convective  transport 

equations  (45) 

=  3.141592653589793 

turbulence  model  constants  for  |l|(=U,V,W,k,e ) 

(1)  transport  quantities  (U,V,W,k,e) 

(2)  y  or  z  in  equation  (42) 
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d,e,  etc.  downstream,  east,  etc.,  control  surfaces 
n,n-l  nth  and  (n-l)th  time  step 


I.  INTRODUCTION 


This  report  is  concerned  with  the  three-dimensional  shear  flow  over  the 
stern  and  in  the  wake  of  a  ship.  A  review  of  the  general  features  of  such 
flows  was  provided  by  Patel  (1982).  On  the  basis  of  available  experimental 
evidence  it  was  observed  that  these  flows  could  not  be  adequately  described  by 
boundary- layer  theory  or  simple  extensions  of  it,  and  recourse  had  to  the  so- 
called  partial ly-parabolic  (or  parabolized)  Navier-Stokes  equations.  The  need 
for  further  experiments  to  provide  the  details  required  for  the  development 
and  verification  of  solution  procedures  was  also  identified.  Over  the  past 
few  years,  much  effort  has  been  devoted  at  many  organizations  to  the  develop¬ 
ment  of  solution  procedures  appropriate  for  the  complex  geometries  of  ship 
hulls,  and  additional  experiments  have  been  conducted  on  a  number  of  different 
ship  and  ship-like  models  to  establish  a  data  base. 

In  the  research  at  The  University  of  Iowa,  a  method  for  the  solution  of 
the  partial ly-parabolic  equations  was  first  developed  and  applied  to  a  variety 
of  trai ling-edge  and  wake  flows,  including  the  flow  over  axisymmetric  bodies 
and  ship-like  three-dimensional  bodies.  This  method  was  described  in  detail 
by  Chen  and  Patel  (1985a)  and  its  applications  have  been  reported  in  several 
publications  [Chen  and  Patel  (1984;  1985a, b),  Patel  and  Chen  (1986a)].  During 
the  course  of  this  development  it  became  clear  that  the  numerical  algorithms 
and  physical  models  incorporated  into  the  method  could  be  employed  equally 
well  in  the  solution  of  the  complete,  fully-elliptic,  Reynolds-averaged  Nav¬ 
ier-Stokes  equations  without  incurring  a  significant  penalty  in  computing 
times  or  storage.  Although,  as  noted  in  Patel  (1982),  the  partial ly-parabolic 
approximations  are  quite  appropriate  for  most  ship  stern  flows,  the  greater 
range  of  applicability  of  the  fully-elliptic  formulation  was  considered  very 
attractive.  For  this  reason,  the  partially-parabolic  method  of  Chen  and  Patel 
(1985a)  has  been  generalized  to  a  fully-elliptic  mode.  The  changes  which  are 
required  in  the  original  formulation  are  described  here.  This  more  general 
method  is  then  applied  to  calculate  the  flow  over  several  ship  forms  for  which 
experimental  data  are  available. 

In  what  follows,  we  shall  first  put  the  present  work  in  perspective  by 
reviewing  the  different  approaches  which  at  being  adopted  to  develop  calcula¬ 
tion  procedures  for  ship  stern  flows.  The  status  of  the  experimental  informa- 
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tion  is  then  reviewed  to  identify  suitable  test  cases  for  the  assessment  of 
calculation  methods.  This  is  followed  by  a  description  of  the  ful ly-el 1 iptic 
method,  and  its  evaluation  by  application  to  the  selected  test  cases.  Detail¬ 
ed  comparisons  are  made  between  the  experiments  and  calculations  with  respect 
to  the  hull  pressure  distribution,  the  mean  velocity  field,  and,  when  avail¬ 
able,  the  turbulence  parameters.  It  is  shown  that  the  numerical  method  is 
capable  of  handling  all  the  geometrical  and  physical  features  of  the  flow. 
Some  improvements  are,  however,  necessary  in  the  modelling  of  the  turbulence, 
particularly  in  the  near-wall  and  near-wake  regions. 

II.  SHIP  STERN  AND  WAKE  FLOW 


I 1.1  DEVELOPMENTS  IN  THEORY 

An  accurate  prediction  of  the  flow  at  the  stern  is  of  great  practical 
interest  in  the  determination  of  ship  resistance,  in  the  design  of  propellers 
and  appendages,  and  in  the  determination  of  the  ensuing  wake.  The  failure  of 
conventional  boundary- layer  calculation  methods,  of  the  type  which  have  been 
used  with  considerable  success  on  aerodynamic  configurations,  to  predict  the 
thickening  of  the  viscous  flow  over  the  stern  was  demonstrated  quite  dramati¬ 
cally  by  the  results  presented  at  the  SSPA-ITTC  Workshop  [Larsson  (1981)]. 
Figure  1,  adapted  from  that  reference,  shows  the  contours  of  the  axial  compo¬ 
nent  of  velocity  at  a  transverse  section  close  to  the  stern  of  two  double 
models  tested  in  wind  tunnels.  There  are  many  reasons  for  the  observed  dis¬ 
agreement  among  the  calculations  with  different  methods  and  experiments. 
Among  these  are:  failure  of  the  boundary- layer  approximations,  inadequacies  of 
the  empirical  inputs  to  effect  closure  of  the  equations,  and,  of  course,  inac¬ 
curacies  in  the  numerical  schemes. 

As  noted  in  the  recent  Report  of  the  Resistance  and  Flow  Committee  of  the 
18th  ITTC  [referred  to  here  as  ITTC  (1987)],  methods  for  the  calculation  of 
stern  flows  have  evolved  in  two  general  directions.  One  of  these  involves 
generalization  and  extension  of  the  thin  boundary- layer  methods,  of  both 
integral  and  differential  types,  to  include  such  factors  as  changes  in  coor¬ 
dinate  metrics  normal  to  the  surface,  normal  variation  of  pressure,  and  inter¬ 
action  between  the  viscous  and  inviscid  flow  regions.  The  development  of 
these  so-called  higher-order  integral  methods  for  thick  boundary  layers  has 
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bejn  pursued  by  Nagamatsu  (1985),  Soejima  (1985),  and  Toda  et  al.  (1985), 
among  others.  These  methods  involve  numerous  approximations  with  respect  to 
the  terms  which  are  retained  in  the  equations,  velocity-profile  assumptions, 
friction  formulas,  and  auxiliary  closure  equations.  Similar  attempts  to 
generalize  differential  thin  boundary-layer  methods  have  been  made  by  Soejima 
(1983)  and  Lee  (1985)  but  they  have  not  reached  the  level  of  development  of 
the  simpler  integral  methods  due  to  the  difficulties  of  coupling  such  calcula¬ 
tions  with  the  external  invicid  flow  to  obtain  matched  and  converged  solu¬ 
tions. 

The  alternative  approach  to  the  calculation  of  stern  flows  involves 
numerical  solutions  of  the  complete  Reynolds-averaged  Navier-Stokes  equations, 
or  the  somewhat  less  general  partial ly-parabolic  equations,  in  domains  which, 
in  principle,  encompass  the  viscous  as  well  as  the  inviscid  flow.  Initial 
development  of  such  methods  was  restricted  to  two-dimensional  and  axisymmetric 
shapes,  but  applications  to  three-dimensional  problems  are  now  beginning  to  be 
made.  Patel  and  Chen  (1986a)  presented  a  review  of  these,  and  other  methods 
of  the  type  discussed  above,  for  the  simpler  case  of  axisymmetric  flow  over 
the  tail  and  in  the  wake  of  bodies  of  revolution.  They  concluded  that  the 
numerical  complexities  of  the  traditional  iterative  approach  to  viscous-invis- 
cid  interactions  is  fast  approaching  those  of  the  global,  large-domain  methods 
which  involve  fewer  assumptions.  However,  comparative  study  of  interactive 
and  noninteractive  procedures  for  the  solution  of  the  partially-parabolic 
viscous-flow  equations  made  by  Stern  et  al.  (1986)  suggests  that  both  approa¬ 
ches  lead  to  satisfactory  results  with  comparable  computing  effort. 

Although  experience  with  global  numerical  methods  in  three-dimensional 
flows,  in  general,  and  for  ship  hulls,  in  particular,  is  still  quite  limited, 
different  approaches  are  being  followed  to  develop  such  methods.  This  is 
evidenced  by  the  recent  work  of  Broberg  and  Larsson  (1984),  Chen  and  Patel 
(1984,  1985a, b),  Hoekstra  and  Raven  ( 1985a, b),  Huang  and  Zhou  (1985),  Ito  and 
Mori  (1985),  Janson  and  Larsson  (1985),  Kodama  (1985,1987),  Raven  and  Hoekstra 
(1985),  Tzabiras  (1984,  1985),  and  Stern  et  al.  (1986).  These  methods  are  at 
various  stages  of  development  and  differ  quite  substantially  from  one  another 
in  many  important  respects.  The  principal  differences  stem  from:  (a)  coordi¬ 
nates  used  and  methods  employed  to  construct  them;  (b)  approximations  intro¬ 
duced  in  the  Reynolds-averaged  Navier-Stokes  equations,  leading  to  the  so- 


called  thin-layer  equations,  or  the  partially-parabolic  equations;  (c)  the 
turbulence  model  employed,  ranging  from  algebraic  eudy-viscosity  to  multi¬ 
equation  models,  and  the  treatment  of  the  boundary  conditions;  (d)  formulas 
used  to  discretize  the  differential  equations  and  the  numerical  methods  used 
to  solve  them;  (e)  in  the  size  of  the  solution  domain  and  coupling  with  the 
inviscid  flow,  if  any;  and  (f)  in  the  manner  in  which  the  solution  is  initi¬ 
ated,  e.g.  using  inviscid  solution,  starting  with  a  boundary-layer  solution, 
marching  in  time  from  rest,  etc.  A  detailed  discussion  of  these  differences 
is  well  beyond  the  scope  of  this  report.  However,  we  shall  address  some  of 
them  in  subsequent  sections  as  we  describe  the  particular  procedures  adopted 
in  the  present  development. 

Most  of  the  methods  mentioned  above  have  been  developed  and  applied  thus 
far  to  the  case  of  double  models.  Thus,  the  free  surface  is  considered  to  be 
flat  and  treated  as  a  plane  of  symmetry.  Extension  of  these  methods  to  treat 
the  effects  of  the  free  surface  have  not  yet  received  much  attention  largely 
because  the  difficulties  are  even  greater  than  those  encountered  within  the 
framework  of  thin  boundary-layer  theory  [see  Stern  (1985,  1986)].  Such  exten¬ 
sions  may  evolve  along  two  different  directions.  In  one,  a  viscous-flow 
method  may  be  combined  with  an  inviscid-f low  method  in  an  interactive  mode  to 
take  advantage  of  the  well  developed  techniques  of  classical  ship-wave  theory. 
The  alternative  approach,  which  does  not  separate  the  viscous  and  wave  ef¬ 
fects,  would  require  the  satisfaction  of  the  proper  boundary  conditions  at  the 
free  surface  which  is  itself  determined  as  a  part  of  the  solution.  This 
latter  approach  has  been  pursued  by,  among  others,  Miyata  et  al.  (1984,  1985, 
1986)  and  Hino  (1987).  In  both  cases,  however,  much  further  work  is  needed  to 
realistically  and  accurately  account  for  the  boundary  layer  and  turbulent  flow 
effects. 

I 1.2  A  GUIDE  TO  EXPERIMENTS 

Measurements  in  ship  boundary  layers  and  in  ship  wakes  have  been  made 
over  many  years.  The  review  of  Patel  (1980)  for  the  Stanford  Conference  on 
Complex  Turbulent  Flows,  which  emphasized  turbulence  modeling,  indicated  that 
none  of  the  data  available  could  be  regarded  as  complete  enough  for  the  vali¬ 
dation  of  stern  and  wake  flow  calculation  methods  because  che  measurements 
were  restricted  to  the  mean  flow.  Among  the  most  detailed  mean-flow  measure- 
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merits  on  ship  hulls  were  those  of  Larsson  (1974)  on  the  SSPA  720  Liner  and 
Hoffman  (1976)  on  the  HSVA  Tanker.  Both  experiments  were  conducted  with 
double  models  in  wind  tunnels.  These  two  data  sets  were  in  fact  used  as  test 
cases  in  the  1980  SSPA-ITTC  Workshop  on  ship  boundary  layers  (Larsson,  1981), 
the  results  of  which  were  mentioned  above  in  connection  with  Figure  1.  These 
experiments  continue  to  be  among  the  most  detailed  available  to  date,  and 
their  value  has  been  greatly  increased  by  the  subsequent  turbulence  measure¬ 
ments  of  Lofdahl  (1982)  [see  also  Lofdahl  and  Larsson  (1984)]  on  the  SSPA 
Liner,  and  the  extensive  wake  measurements  of  Wieghardt  and  Kux  (1980)  and 
Wieghardt  (1982,  1983),  and  turbulence  measurements  of  Knaack  (1984)  and 
Knaack,  Kux  and  Wieghardt  (1985)  in  the  wake  of  the  HSVA  tanker.  The  data  on 
these  two  hull  forms,  therefore,  are  detailed  enough  to  test  several  critical 
aspects  of  ship  stern  and  wake  flow  calculation  methods. 

The  Cooperative  Experimental  Program  of  the  ITTC  Resistance  and  Flow 
Committee  provided  an  impetus  to  the  establishment  of  a  comprehensive  data 
base  on  the  flow  around  ship  hulls  against  which  emerging  theoretical  methods 
could  be  tested.  While  the  focus  of  this  program  was  on  quantities  of  princi¬ 
pal  concern  in  tankery,  it  also  led  to  some  detail  measurements  in  the  viscous 
flow.  The  present  status  of  the  program  is  reviewed  in  ITTC  (1987).  Of  the 
four  hulls  which  were  initially  selected  for  the  program,  which  included  the 
HSVA  Tanker  mentioned  above,  measurements  in  the  boundary  layers  and  wakes  of 
two  have  been  reported  to  date.  These  are  the  Wigley  parabolic  form  and  the 
Series  60,  Cg  =0.60  form.  In  spite  of  the  concerted  international  effort  over 
a  period  of  years,  not  all  of  the  available  data  are  suitable  as  test  cases 
for  numerical  methods  of  the  type  described  in  the  previous  section.  "s 
pointed  out  in  ITTC  (1987),  many  experiments  were  confined  to  a  few  stations 
over  the  stern  and  in  the  very  near  wake  and,  therefore,  it  is  difficult  to 
establish  proper  initial  and  boundary  conditions  for  them.  Some  involve  large 
scatter  and  uncertainty  due  to  blockage  and  model  attitude.  Also,  quite 
different  measurement  locations  and  coordinates  have  been  employed  by  differ¬ 
ent  experimenters  even  for  tests  on  the  same  hull  form.  While  this  precludes 
direct  comparisons  among  data  sets  from  different  experiments  and  test  facil¬ 
ities,  some  of  the  data  sets  are  sufficiently  well  documented  for  use  iri 
testing  the  capabilities  of  calculation  methods  to  predict  certain  important 
aspects  of  the  flow. 
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In  addition  to  the  measurements  on  the  Wigley  and  Series  60  hulls  made 
under  the  auspices  of  the  ITTC  program,  data  have  been  obtained  in  the  stern 
flow  of  other  ship  hulls  and  ship-like  forms.  Fukuda  and  Fujii  (1985)  have 
reported  measurements  on  three  hulls,  including  a  Series  60,  Cg=0.8  form  and 
an  elongated  model  derived  from  it  by  adding  a  parallel  middle  body  to  study 
scale  effects,  by  Hotta  and  Hatano  (1985)  on  a  tanker  model,  and  by  Hinatsu 
and  Takeshi  (1985)  on  two  other  hulls.  Of  these,  the  first  were  conducted  on 
double  models  in  a  wind  tunnel,  the  second  used  a  free-surface  model  in  a 
water  channel,  and  the  last  was  performed  in  a  towing  tank.  Only  the  first 
two  included  turbulence  measurements.  Among  the  experiments  on  ship-like 
bodies  are  those  of  Huang  et  al.  (1983)  on  two  models  of  elliptic  cross  sec¬ 
tion  with  the  same  sectional  area  distributions  as  that  of  one  of  the  axisym- 
metric  bodies  which  were  tested  earlier. 

Although  the  foregoing  summary  of  experiments  gives  the  impression  that 
there  is  now  a  substantial  body  of  experimental  information  on  stern  flows 
over  a  variety  of  shapes,  in  reality  the  available  data  are  not  sufficiently 
complete  or  extensive  enough  to  be  used  as  comprehensive  test  cases  for  calcu¬ 
lation  methods.  In  fact,  this  deficiency  is  not  restricted  to  ship  flows.  It 
is  shared  by  practically  all  measurements  in  three-dimensional  flows  because 
such  measurements  are  time  consuming  and  expensive,  and  the  experimental 
techniques  themselves  are  not  well  developed.  Also,  the  experiments  are 
rarely  conducted  for  the  purpose  of  validating  all  aspects  of  any  particular 
calculation  method.  The  choice  of  test  cases  for  this  purpose  is,  therefore, 
a  difficult  one.  We  shall  address  this  issue  in  a  later  section. 

I 1.3  SCOPE  OF  THE  PRESENT  WORK 

The  purpose  of  this  report  is  two-fold:  (a)  to  describe  the  generaliza¬ 
tion  of  the  partially-parabolic  method  of  Chen  and  Patel  (1985a)  to  a  fully- 
elliptic  capability,  and  (b)  to  assess  the  performance  of  the  method  by  com¬ 
parisons  between  calculations  and  experiments.  The  first  aspect  is  relatively 
straightforward  because  the  partially-parabolic  method  was  described  in  con¬ 
siderable  detail  in  the  above  reference,  and  the  changes  in  the  numerical 
structure  of  the  method  are  not  extensive.  In  the  interest  of  completeness, 
however,  the  essential  parts  of  the  method  are  reviewed.  The  method  is  de¬ 
scribed  in  the  next  section  (Section  III). 
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The  second  task  is  more  formidable  for  several  reasons.  First,  it  is 
important  to  identify  the  most  critical  and  the  most  successful  aspects  of  the 
overall  method  by  evaluating,  to  the  extent  possible,  the  performance  of  each 
of  the  many  components  of  the  method.  This  involves  the  selection  of  simple 
test  cases  in  which  a  particular  feature  of  the  method  can  be  evaluated. 
Fortunately,  the  present  method  has  been  subjected  to  such  tests  and  the 
results  have  been  reported  in  the  literature.  For  example,  the  performance  of 
the  k-e  turbulence  model,  together  with  the  two-point  wall  functions  approach, 
features  which  are  retained  in  the  present  method,  has  been  examined  in  the 
previous  applications  of  the  partially-parabolic  method,  and  the  capability  of 
the  elliptic  numerical  scheme  to  handle  separation  and  reattachment  in  axisym- 
metric  laminar  and  turbulent  flows  was  demonstrated  by  Patel  and  Chen  (1986b) 
and  Chen  and  Patel  (1987b,  1988),  respectively.  Secondly,  it  is  necessary  to 
evaluate  the  complete  methodology  by  comparisons  with  experimental  data  on  a 
wide  variety  of  realistic  ship  forms.  Here,  the  previous  experience  is  some¬ 
what  limited.  The  partially-parabolic  method  was  first  compared  with  data  on 
the  relatively  simple  ship-like  bodies  of  elliptic  cross  section  in  Chen  and 
Patel  (1984,  1985a)  and  with  experiments  on  the  Wigley  and  SSPA  hulls  in  Chen 
and  Patel  (1985b).  Although  satisfactory  agreement  was  observed  in  all  cases, 
some  difficulties  were  noted,  particularly  in  the  last  reference,  with  regard 
to  grid  refinement  and  turbulence  modeling.  That  reference  also  pointed  out 
yet  another,  and  equally  important  difficulty  of  carrying  out  detailed  compar¬ 
isons  between  calculations  and  experiment.  This  relates  to  the  different 
coordinates  employed  in  the  two,  and  the  loss  of  accuracy  resulting  from  the 
required  interpolation.  While  these  difficulties  cannot  be  readily  resolved, 
they  emphasize  the  need  to  make  comparisons  with  a  number  of  data  sets,  col¬ 
lected  in  different  ways  and  on  different  models,  to  provide  a  more  complete 
account  of  the  capabilities  of  the  numerical  method.  Special  attention  is 
therefore  paid  in  the  present  work  to  examine  the  available  data  summarized  in 
the  previous  section  to  identify  potential  test  cases.  The  rational  for 
selecting  five  hull  forms,  namely  the  Wigley  parabolic  ship,  the  SSPA  720 
Cargo  Liner,  the  HSVA  Tanker,  the  SR107  Ore  Carrier,  and  the  Series  60,  C^  = 
0.6  form,  will  be  discussed  along  with  the  corresponding  calculations  and 
comparisons  with  data  in  Section  IV. 
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III.  CALCULATION  STRATEGY  AND  METHOD 


III.l.  EQUATIONS  AND  COORDINATES 

For  three-dimensional  flows  involving  complex  geometries,  it  is  desirable 
to  employ  body-fitted  coordinate  systems  so  that  the  flow  in  the  wall  layer 
can  be  accurately  resolved  with  a  reasonable  number  of  grid  points.  Once  such 
a  coordinate  system  is  selected  for  a  given  geometry,  there  remains  the  task 
of  formulating  the  equations  of  motion  in  that  system.  Two  different  approa¬ 
ches  can  be  adopted  for  this  purpose.  One  of  these  uses  what  may  be  termed 
"partial  transformations",  in  which  only  the  independent  coordinate  variables 
are  transformed,  leaving  the  dependent  variables  (i.e.  velocity  components)  in 
a  preselected  coordinate  system  in  the  physical  domain.  This  approach,  which 
has  been  used  by  Chen  and  Patel  ( 1985a, b),  among  others,  has  the  advantages 
that  the  resulting  equations  are  relatively  simple  and  the  results  can  be 
readily  interpreted.  Since  the  velocity  vectors,  in  general,  do  not  align 
with  the  coordinate  directions,  this  approach  may  lead  to  increased  numerical 
diffusion  when  the  angles  between  the  velocity  components  and  coordinate  sur¬ 
faces  become  large.  The  alternative  is  to  transform  the  equations  completely, 
including  the  independent  as  well  as  the  dependent  variables.  This  approach 
has  been  used  by,  among  others,  Richmond  et  al.  (1986),  Stern  et  al.  (1986), 
and  Ogawa  and  Ishiguro  (1987).  The  use  of  contravariant  velocity  components 
in  such  a  complete  transformation  allows  a  much  more  accurate  resolution  of 
the  boundary-layer  flow  near  a  solid  surface.  However,  the  fully-transformed 
equations  involve  many  geometric  coefficients  and  their  higher-order  deriva¬ 
tives.  This  not  only  leads  to  increased  computer  storage  requirements  but 
also  can  adversely  affect  the  solution  of  the  flow  equations  if  the  coeffici¬ 
ents  are  not  evaluated  accurately.  In  many  practical  applications,  it  is  not 
necessary  to  use  the  complete  transformations  if  the  basic  coordinate  systems 
are  chosen  carefully  so  as  to  avoid  large  skew  angles  between  velocity  compon¬ 
ents  and  the  faces  of  the  computational  cell. 

In  the  present  study,  which  is  concerned  with  external  flows  past  ship 
forms,  the  restrictions  associated  with  the  partial  transformations  can  be 
easily  alleviated  by  choosing  cylindrical  polar  coordinates  as  the  basic 
coordinate  system  in  the  physical  plane.  A  detailed  description  of  such 
partial  transformations  and  the  corresponding  equations  was  given  in  Chen  and 
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Patel  (1985a).  For  the  sake  of  completeness,  we  shall  briefly  describe  these 
in  the  following. 

Consider  the  equations  of  motion  in  cylindrical  coordinates  (x,r,0)  for 
unsteady,  three-dimensional,  incompressible  flow.  The  exact  Reynolds-averaged 
equations  of  continuity  and  momentum  of  the  mean  flow,  in  dimensionless  form, 
are, 
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where  x  =  jj,  r  =  j-,  0  are  the  dimensionless  coordinates  normalized  by  a  char¬ 
acteristic  length  L,  and  t  is  the  time  normalized  by  L/UQ.  U,  V,  W  are, 
respectively,  the  longitudinal,  radial  and  circumferential  components  of  mean 
velocity  normalized  by  the  characteristic  velocity  U0.  P  is  the  pressure 
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normalized  by  pUq  .  Re  =  UQL/v  is  the  Reynolds  number  defined  in  terms  of  U0> 

L  and  molecular  kinematic  viscosity  v.  The  barred  quantities  uu,  uv,  etc., 
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are  the  Reynolds  stresses,  normalized  by  U  . 

In  the  present  study,  the  two-equation  k-e  turbulence  model  is  used  to 
model  the  Reynolds  stresses.  Each  stress  is  related  to  the  corresponding  mean 
rate  of  strain  by  an  isotropic  eddy  viscosity  v  as  follows: 
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The  eddy  viscosity  is  related  to  the  dimensionless  turbulent  kinetic  energy  k, 
and  its  rate  of  dissipation  e,  by 


where  c^  is  a  constant,  and  k  and  e  are  governed  by  the  convective  transport 
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where  G  is  the  turbulence  generation  term: 
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The  effective  Reynolds  number,  R  ,  is  defined  as 
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Where  <t>  =  (U,V,W,k,e).  The  constants  in  these  equations  were  taken  as 

Cu  =  0.09,  Cel  =  1.44,  Ce2  =  1.92,  °u  =  orv  =  °w  =  °k  =  h0,  °e  =  1*3*  We  note 

again  that  v^,  k,  and  e  are  all  made  dimensionless  by  the  characteristic 

velocity  and  length  scales,  U0  and  L,  respectively. 

Using  equations  (5,  6  and  10),  the  momentum  equations  (2)  through  (4)  and 
the  turbulence-model  equations  (7)  and  (8)  can  be  written 
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a,  r  30  1  r  36 
k 


1  2 

F" 7  k  -  G  +  e  =  0 
Kk 


3v 


3v, 


3v , 


3£./„  i_  L_£\  .  /»  J _ ttle./u  _L  _i\  1 

3T  +  a  3x  ’  3x  +  '*  a  3r  7  3r  +  a  r  38  '  r 


3e 

30 


J—  v2e  -  f  —  R  +  r  —  =n 
Re  Lel  k  b  +  'e2  k  U 


(12) 


(13) 


(14) 


(15) 


It  is  convenient  to  rewrite  these  transport  equa  ions  (11)  through  (15)  in  the 
following  general  form 

V2*  =  R*  [<U  -  ytiX)  ♦„♦<*-  Vt>r!  *r  ♦ 

(H  -  <*♦  r  vt.o>  <F*o>  *  +  *♦  1161 

where  9  again  represents  any  one  of  the  convective  transport  quantities:  U,  V, 

W,  k  or  e,  and  the  subscripts  x,r,0  denote  derivatives.  The  corresponding 

coefficients  b,,  c.  and  d.  are 
9  9  9 
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bu  =  2> 


CU  =  1* 


du  =  1 


bv  =  1. 

V’ 

bk  =  o“* 


be  =  a  * 
e 


cy  =  2, 

CW  =  lf 
c  =J- 

k  V 


'e  a  * 
e 


dv  =  i 

d«“2 
d  -J- 

k '  V 


e  a  » 
e 


(17) 


and  the  source  functions  s^  for  U,  V,  W,  k  and  e  are,  respectively, 
su  *  RU  tPx  +3  V  vt,rV  <r  vt,9> 


(18a) 


2 

s v  =  Ry  [pr  +  ^  kr"  vt  xUr  "  ^r  vt  0^Wr  “  r^  “  r~^ 


2  V 

Twe  +T 


(18b) 


SW  =  RW  '-r  p0+  3  r  k0  +  r  '  vt ,x*r  V 


/ 1  w  W« 
-  vt,r  r  V  r* 


,1  *  ,2V , ,  2_  W 

"  ^r  vt,0 ^  ^  “  2  V0  +  2 

r  r 


(18c) 


Sk  = 


se  = 


Rk  (G  -  e) 


Re  (Cel  k  G 


Ce2  k  ) 


(18d) 

(18e) 


Rations  (11)  through  (15)  are  coupled,  nonlinear,  partial  differential 
nations  and,  together  with  the  continuity  equation  (1),  are  sufficient,  in 
Mnciple,  to  solve  for  the  six  unknowns  p,  U,  V,  W,  k  and  e  when  proper 
litial  and  boundary  conditions  are  specified. 
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(a)  Body-Fitted  Coordinates 

Three-dimensional  geometries,  such  as  ship  hulls  are  usually  quite  com¬ 
plex  and  cannot  be  conveniently  described  by  simple  orthogonal  coordinate 
systems.  It  is,  therefore,  desirable  to  introduce  analytic  or  numerical 
coordinate  transformations  which  simplify  the  computational  domain  in  the 
transformed  plane  and  facilitate  applications  of  the  boundary  conditions.  In 
the  present  study,  we  shall  use  a  numerically-generated,  body-fitted  coordin¬ 
ate  system,  since  it  offers  the  advantages  of  generality  and  flexibility  and, 
most  importantly,  transforms  the  computational  domain  into  a  simple  rectangu¬ 
lar  region  with  equal  grid  spacing. 

In  the  numerical  grid-generation  technique,  we  seek  a  coordinate  system 

for  the  numerical  analysis  of  the  flow  in  the  domain  D  shown  in  Figure  2. 

This  domain  is  bounded  by  an  arbitrary  hull  surface  S,  the  ship  centerplane  C, 

the  free  surface  or  water  plane  W,  the  upstream  and  downstream  sections  A  and 

B,  respectively,  and  an  external  boundary  £.  Section  A  may  be  located  at  a 

hull  section  where  the  boundary  layer  is  thin  to  avoid  calculation  of  the  flow 

over  the  bow  or  it  may  be  placed  far  ahead  of  the  ship  if  the  bow  is  to  be 

included.  The  downstream  boundary  B  may  be  placed  at  a  section  in  the  far 

wake.  The  choice  of  the  external  boundary  £  is  also  arbitrary.  It  could  be 

far  away  from  the  hull  or  it  could  coincide  with  the  walls  of  a  towing  tank  or 

wind  tunnel.  The  basic  idea  of  a  boundary-conforming  curvilinear  system  is  to 

find  a  transformation  such  that  the  boundary  surfaces  of  the  physical  domain  D 

in  cylindrical  or  in  any  other  basic  orthogonal  coordinate  system,  say 
1  2  3 

(x  ,  x  ,  x  ),  are  transformed  into  boundaries  of  a  simple  rectangular  domain 
in  the  computational  space  (5,  n,  t)  shown  in  Figure  3. 

With  the  values  of  the  curvilinear  coordinates  specified  on  the  boundar¬ 
ies  of  0,  it  remains  to  generate  the  values  of  these  coordinates  in  the  inter¬ 
ior  of  D.  This  is  a  boundary-value  problem  in  the  physical  field  with  the 
curvilinear  coordinates  (£,  n,  c)  as  dependent  variables  and  the  orthogonal 
coordinates  (x  ,x  ,x  )  as  the  independent  variables,  with  boundary  conditions 
specified  on  the  curved  boundaries.  Thus,  a  system  of  elliptic  partial  dif¬ 
ferential  equations  can  be  used  to  generate  the  coordinates  since  the  field 
solution  of  such  a  system  is  determined  entirely  by  the  boundary  conditions. 
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However,  the  elliptic  system  must  be  chosen  such  that  it  precludes  the  occur¬ 
rence  of  extrema  in  the  interior  of  the  domain  and  assures  a  one-to-one  map¬ 
ping  between  the  physical  and  transformed  planes. 

For  the  general  three-dimensional  but  simply-connected  domain  of  interest 
here,  a  set  of  Poisson  equations  of  the  form 


v2?  =  f 1  (? ,  n,  O 
v2n  =  fZU,  n ,  O 
v2?  =  f3 (5 ,  n,  O 
or  simply 

2  i  i 

vV  =  f  ,  i  =  1,  2,  3 

1  2  3 

with  5  =  5,  5  =  n,  5  =  c 


(19) 


(20) 


2 

may  be  taken  as  the  coordinate  generating  system.  Here,  v  is  the  Laplacian 
operator  in  orthogonal  coordinates  x1.  The  nonhomogeneous  source  functions  f1 
may  be  assigned  appropriate  values  to  yield  the  desired  concentration  of 
coordinate  surfaces.  The  choice  of  these  functions  for  specific  applications 
will  be  discussed  later  in  Section  111.2(a).  Equation  (20)  is  subject  to 
either  Dirichlet  or  Neumann  boundary  conditions  on  the  boundary  surfaces, 
which  are  surfaces  of  constant  S1. 

Since  it  is  desirable  to  perform  all  numerical  computations  in  the  trans¬ 
formed  (S,  n,  c)  plane  with  equal  grid  spacing,  i.e..  a?  =  An  =  a?  =  1, 

equation  (20)  is  cumbersome  to  use.  It  is  more  convenient  to  invert  it  and 
solve  for  the  orthogonal  coordinates.  In  other  words,  the  dependent  and 

independent  variables  are  interchanged  so  that  the  orthogonal  coordinates 
12  3 

(x  ,x  ,x  )  in  the  physical  plane  become  the  dependent  variables,  with  the 
curvilinear  coordinates  (5,  n,  s)  as  the  independent  variables.  The  boundary- 
value  problem  in  the  transformed  field  then  involves  generating  the  values  of 
the  orthogonal  coordinates  x1  =  x1  (?,  n,  c)  in  the  interior  from  the  speci¬ 
fied  boundary  values  of  x1  on  the  rectangular  boundary  surfaces  of  the  trans¬ 

formed  field.  Since  the  boundaries  in  the  transformed  plane  are  all  rectangu- 
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lar  (constant  S,  n,  or  s  plane),  these  computations  are  carried  out  on  a  cubic 
grid  regardless  of  the  shape  of  the  physical  boundaries. 

In  obtaining  the  inverse  transformation  of  equation  (20),  several  general 

12  3 

relations  between  the  physical  (x  ,x  ,x  )  and  the  transformed  (5,  n,  s)  coor¬ 
dinates  are  required.  The  basic  expressions  may  be  found  in  many  reference 
books,  for  example  Aris  (1971)  and  Lass  (1975),  although  some  of  the  relations 
are  not  explicitly  given.  A  summary  of  some  important  relations  are  given  in 

Appendix  I  of  Chen  and  Patel  (1985a)  with  specific  reference  to  the  transfor- 

12  3 

mat  ion  between  an  orthogonal  coordinate  system  (x  ,x  ,x  )  and  general  coord m- 
12  3 

ates  (S  ,  S  ,  S  )  =  (Sr  n,  ?).  The  latter  are  not  necessarily  orthogonal. 
With  these  transformation  formulae,  equations  (20)  become: 


vV  - 


hlh2h3  ax1 


h  h 
(-2-2.1 
1  h  ) 


v2x2  - 


hlh2h3  3x2  h2 


1 


(21) 


72x3 


_ ] _ _3_ 

hlh2h3  3x3 


where 


11  a 

g  — o  +  g 

as 


+  2  g 


23  a' 


2 

22  a^ 

2 

x  33  _il 

1  2 

+  9  2  1 

3n 

as 

_1 

a  ,2  a 

S  +  f 

as  +  f  an 

2  2 
9  12  _iL_  .  9  13  _L_ 
29  3?3n  +  29  3^3C 


.3  a 


as 


(21a) 


h-j  are  the  metric  coefficients  in  the  chosen  orthogonal  coordinates  x1,  and 
g1J  is  the  conjugate  metric  tensor  in  the  transformed  coordinates  s'(=S,n,s). 
Note  that,  for  cylindrical  polar  coordinates  (x,  r,  0),  h-j  are  (1,  1,  r), 
respectively. 

Equations  (21)  can  be  solved  numerically  in  the  transformed  domain 
(S,n,s)  when  proper  boundary  conditions  are  specified  on  all  boundary 
surfaces  (i.e.,  constant  S,  n  and  s).  If  f  =  f  =  f  =  0,  the  transformation 
is  said  to  be  homeomorphic.  In  general,  however,  non-zero  values  are  assigned 
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to  these  functions  to  exercise  control  over  the  grid  distribution.  Solutions 
of  equations  (21)  to  obtain  numerical  grids  for  particular  shapes  are  present¬ 
ed  in  Section  111.2(a). 


(b)  Reynolds  Equations  in  Transformed  Coordinates 

Although,  as  noted  earlier,  it  is  possible  to  transform  both  the  indepen¬ 
dent  (x\  x^,  x^,  t)  and  the  dependent  (U,  V,  W,  k,  e)  variables  to 
the  (5,n,s)  coordinates,  we  will  consider  only  the  transformation  of  the  inde¬ 
pendent  variables,  leaving  the  velocity  components  U,  V,  W  in  the  original 
(x\  x2,  x^)  coordinates  in  the  physical  plane.  Also,  as  mentioned  earlier, 
the  cylindrical  polar  coordinates  (x,  r,  0)  appear  to  be  most  convenient  for 
the  description  of  the  flow  field  around  practical  ship  forms.  Thus,  in  this 
study,  equations  (1)  and  (16)  [i.e.,  (11)  through  (15)]  will  be  used  as  the 
basic  equations  in  the  physical  plane  to  derive  the  equations  Ip  the  trans¬ 
formed  domain  (C,  n,  ?).  Transformed  equations  resulting  from  other  ortho¬ 
gonal  coordinates,  which  may  be  useful  in  other  applications,  are  given  in 
Appendix  I  of  Chen  and  Patel  (1985a).  Using  these  general  transformation 
relations,  equation  (16)  for  an  unsteady  three-dimensional  flow  can  be  written 
in  the  following  form: 


11  2?  33  12  13  23 

9  V  9  ♦nn+  *  2g  W  *9  ♦«+  2g 


V  f2v  fV  f  ^  <»!v  biV  bk> 


+  f 


+  %2  <b2  »c  +  b2  V  b2  V  +  %  <b3  V  b3  V  bk)] 


*  Vt  *  S* 


(22) 


where 


ex]  =  U  -  -A  (bl  v,  +  b?  v.  +  b?  v  , 

<J>  J  l  t,c  i  t,n  I  t,< 

*1  -  V  -  ^  (b2  V  b2  Vt,n  +  ^  -,J 

a*  =  w  -  'f  (b3  vt,5+  b3  vt,n  +  b3  vt,t} 


(23) 
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and  b. ,  c,  and  d,  are  as  defined  in  equation  (17).  Hereafter,  the  subscripts 
4>  «p  <f> 

(?,n,t)  on  4>  (=  U,V,W,k,e)  and  v  denote  derivatives. 

Equation  (22)  is  identical  with  equation  (A-84)  in  Appendix  I  of  Chen  and 
Patel  (1985a).  It  can  be  rearranged  into  a  general  convective-transport 
equation  of  the  form 

n  ??  qq 

g  *55+9  *„„+  g  V’ 2\V  2B»V  2C»V  Vt  +  s*  (24> 

where 


1  /h3J  ,  .3  2  .  h3  3v  f3 
-  j  +  b2%  +  b3a^  ~  f 


(25a) 


2B*  =  f  (bla*  *  b¥l  *  b¥¥ 


(25b) 


2C<J>  J  (bl  4>  +  b2  <t>  +  b3a<f> 


(25c) 


V  S*  -  2<An  +  Ac  *  A,) 


(?5d) 


and  the  source  functions  s.  are  given  by 

<p 

5U  *  Ru  [J  (bl  h*  bk+  bl  P;>  +  3j  (blV  blkn+  blV 


'  ^  (b2  vt,{  +  b2  vt,n+  b2  vt,?)(blV  blV  blV 

J 

'  ^  lb3  ut,«+  b3  ”t,n+  b3  WE+  bllJn  +  blV] 


(26a) 


SI  ‘  V"  r  +  J  (b2?5*  b2pr,+  bk’  +  3J  <b2kE+  b2kn+  b2kc  * 


7  <b!“t,{  +  bl  vt,n+  bl  vt,E^b2UE+  blv  b2Uc> 


■  J  (b3  vt,C+  b3  vt,n+  b3  (b2W?+  b2Wn 

+  b|  W? )  -  7]}  +  (bjw?+  b3^+  b|w? )  +  \ 

r 


(26b) 


SW  =  RW  {  r  +  J  ^b3pC  +  b3pn4  b3pC^  +  3J  *b2kS  +  b3kn+  b3kC^ 

‘  ~j2  (blvt,S  +  b1vt,n+  blvt,?)(b3US+  b3Un+  b3*V 

“  J  (b2Vt,5+  b2vt,n+  b2vt,?)!;j  (bX+  b3Vn+  b3V^  '  r1 

2V  .1  .2  ,3  4,  2  ,,  1,.  ,2,.  ,3  .  w 

-  71  (b3vt,c+  b3vt,n+  b3vt^)}  -  7j  (b3V?  +  b3Vn  +  bsV  +  ^ 


(26c) 


Sk  =  "  Rk(G  -  e) 


(26d) 


Se  ’  -  Re  (C£1  t  G  *  C£2  k  > 


(26e) 


G  =  «t  {\  <b]v  biun+  bX)2+  h  (b2V  b2Vn+  b2VC)2 
J  J 

*  2  [J  (b’w5+  b2W„+  b^)  ♦  i]2 

+  ~  <b]v?+  b2Vn  +  bGV?+  b,Ju5  +  b|un+  b^)2 

♦  \  (b]«?+  b2Wnt  bjUjt  b2U„  ♦  b^2 

J 

+  [-j  (bjw5+  b2Wn  +  b^+  b^V^  b2Vn+  b^)  -  ~]2} 


(26f) 


The  radius  r,  the  geometric  coefficients  b!?  and  g1J,  and  the  Jacobian  J 
which  appear  in  the  above  equations  are  functions  of  the  coordinates  only. 
When  either  analytic  or  numerical  transformations  are  employed  to  generate  the 
grid  distribution,  r  =  r(5,  n,  t)  is  known  in  the  transformed  plane.  The 
geometric  coefficients  b^  and  g1J  for  the  present  transformation  are 
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(b  h  = 


'  1  K1  J] 
b1  b2  b3 

.2  h2  K2 

b1  b2  b3 

h3  h3  K3 

bl  b2  b3 


r(rnV  r?en)  rU^-  x^) 
r(rWrcV  r(x?e?-  x^)  x?r?-  x?r? 

r(rcV  rnV  ^nVVr1  Vn'Vs 


99  =  922  933  "  9232 

9922  =  gn  933  -  3132 

gg33  =  gn  §22  -  9i22 

gg1 2  =  9g21  =  gi3  923  '  912933 

gg13  =  gg31  =  932  923  ~  g i 3^22 

gg23  =  gg32  =  g12  g]3  -  g23g11 


where 


2  2  2  2 

911  =  XC  +  r?  +  r  0C 

2  2  2  2 

922  =  %  *  rn  +  r  K 

2  2  2.2 

933  =  xc  +  r5  +  r  6? 


912  ‘  92)  •  *c*n  *  r?rn  +  r  9^ 

913  =  93l  =  Vc  *  Vi  *  rV? 

923  s  932  '  V{  ♦  rnrs  *  r29ne? 


9  =  9n  Sjj  923  +  29l2  9i3  933  -  l933)2  9n  "  (9l3)2  S32 
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r\ 

~  ^912'  933 


(30) 


is  the  determinant  of  the  metric  tensor  g^.  Also,  the  Jacobian  can  be  ex¬ 
pressed  as 


(31) 


It  should  be  recalled  here  that  (U,  V,  W)  represent  the  velocity  compon¬ 
ents  in  the  physical  cylindrical  polar  coordinates  (x,  r,  0).  Equations  (24) 
through  (26),  together  with  the  equation  of  continuity  (1),  which  transforms 
to: 


(b]u  +  b^V  +  b]w)?  +  (b'ju  +  b^V  +  b^  +  (b^J  +  b^V  +  b|w)^=  0 


are  the  Reynolds-averaged  Navier-Stokes  equations  for  unsteady,  three-dimen¬ 
sional  flows. 


I I I. 2.  NUMERICAL  SOLUTION  PROCEDURES 

In  this  section  we  describe  the  numerical  techniques  used  for  the  solu¬ 
tion  of  the  grid-generation  equations  (21),  the  ful ly-el 1 iptic  convective- 
diffusion  equations  (24),  and  the  equation  of  continuity  (32).  Details  of  the 
numerical  method  used  in  grid  generation,  together  with  some  examples,  are 
presented  in  Section  (a).  The  finite-analytic  numerical  method  of  Chen  and 
Chen  (1982,  1984)  and  Chen  and  Patel  (1985a)  is  then  revised  and  extended  to 
solve  the  five  transport  equations  for  <)>  =  U,V,W,k  and  e  with  a  guessed  pres¬ 
sure  field.  This  is  described  in  Section  (b).  The  continuity  equation  is 
used  to  obtain  equations  which  enable  the  determination  of  the  pressure  and 
pressure-correction  fields  from  the  velocity  field.  This  procedure  is  des¬ 
cribed  in  Section  (c).  The  treatment  of  the  boundary  conditions  is  described 
in  Section  (d),  and  the  overall  numerical  solution  algorithm  is  summarized  in 
Section  (e). 
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(a)  Grid  Generation 

The  numerical  curvilinear  coordinate  system  is  constructed  by  solving 
equations  (21).  Because  the  orthogonal  coordinates  x1  in  the  physical  plane 
can  be  selected  arbitrarily  and  the  control  functions  f1  are  independent  of 
these  coordinate,  it  is  possible  to  generate  the  numerical  coordinates  in 
terms  of  the  simplest  Cartesian  coordinates  (x,y,z),  and  then  transform  them 
into  other  orthogonal  coordinates  which  are  chosen  to  specify  the  velocity 
components  in  the  equations  of  motion.  As  noted  earlier,  it  is  desirable  to 
use  cylindrical  polar  coordinates  as  the  basic  coordinate  system  for  the 
present  applications  to  three-dimensional  ship  forms.  It  is  also  desirable  to 
choose  the  transverse  sections  of  the  hull  as  the  constant-C  stations, 
i.e.,  5  =  5(x),  so  that  the  computing  effort  required  for  grid  generation  and 
flow  calculations  can  be  significantly  reduced.  With  this  choice,  equations 
(21)  reduce  to 


g]1x5?+  f]x?  =  0 


911r«+  92%n+  933rV  2912^n+  2g13ryr+  2g23r, 


5n 


ne 


+  fV  +  f2r  +  f^r  -  — 

+  r  r^+  t  rn+  r  r^  r 


11  22  33  1?  13  23 

9  V  9  °nn+  9-%+  2g‘%n+  2gl\?+  2g236 


(33) 


nc 


1  2  3 

+  f  e  +  t  8  +  f  e  =  o 
€  n  C 


relating  the  numerical  coordinates  (?,n,c)  to  the  cylindrical  coordinates 
(x,r,9)  in  the  physical  plane.  Equivalently,  these  equations  can  also  be 
written  in  Cartesian  coordinates  (x,y,z)  as 

g]1x5?+  f]x^=  0 

11  22  33  12  13  23 

9  -%+  9  ynn+  9  2g  y?n+  2g  y??+  2g  ynC 

♦  f,y?  +  fV  fV  o  ,34, 
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Jl,  X  n22,  X  n33,  x  oJ2,  x  9n137  x  On23, 

9  z5!+  9  znn+  9  zcc4  2g  z?n+  29  ZU+  29  zn? 

+  f]z,+  f2z  +  f3z  =  0 
Cut 

where  y  =  r  sin6,  z  =  r  cose,  and  6  is  the  circumferential  or  girthwise  angle 
measured  from  the  keel  to  the  waterline.  The  control  functions  f1  are  the 
same  in  both  coordinate  systems,  and,  in  principle,  equations  (33)  and  (34) 
yield  the  same  numerical  coordinates  if  the  same  control  functions  are  em¬ 
ployed.  Numerically,  however,  the  xn-coordinates  used  to  specify  the  hull 
geometry  influence  the  accuracy  of  the  calculated  numerical  coordinates  due  to 
truncation  errors.  Since,  for  ship  hulls  the  variation  of  the  surface  coor¬ 
dinates  (y,z)  is  much  smoother  than  that  of  (r,e)  it  is  preferable  to  use 
equations  (34)  instead  of  (33).  The  numerical  coordinates  thus  obtained  are 
then  transformed  to  (x,r,e)  which  are  used  to  specify  the  velocity  components 
in  the  equations  of  motion. 

For  the  numerical  solutions,  it  is  convenient  to  rewrite  equations  (34) 
in  the  form 

g]1(x55-  2ax^ )  =  0 

g]1(y?r  2ay5)  +  g22(ynn-  2byn)  +  g33(y^-  2cy^) 

+  2g12y5n+  2g13ycc  +  2g2V  0 

g11(z?r  2az5)  +  g22(znn-  2bZri)  +  g33(z?r  2cz?) 

1?  13  23 

+  2g  Sn +  2g,%,  +  2g2%c=  0 

where  a,  b  and  c  are  modified  control  functions  defined  by 


(35a) 

(35b) 

(35c) 
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2b  =  - 


(36) 


2c  =  - 


2g 

f 


2g 


2 

22 

3 
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In  order  to  solve  equations  (35),  it  is  necessary  to  prescribe  the  boundary 
conditions  and  the  control  functions.  The  boundary  conditions  are  determined 
by  the  preselected  size  of  the  solution  domain  and  the  desired  number  of  grid 
points  in  the  axial  (5),  radial  (n)  and  circumferential  (c )  directions. 


In  all  the  calculations  considered  here,  the  first  station  5  =  1,  is 
located  at  x  =  0.3  where  the  boundary  layer  is  thin,  and  the  last  station  £  = 
LL  is  placed  in  the  far  wake  where  upstream  viscous  diffusion  is  negligible. 
In  the  radia1  direction,  there  are  MM  points  with  n  =  1  corresponding  to  the 
hull  surface  or  wake  centerplane,  and  n  =  MM  being  the  exterior  boundary.  Tne 
latter  is  placed  typically  at  a  distance  of  the  order  of  one  ship  length  from 
the  hull  for  calculations  corresponding  to  an  unrestricted  stream. 
Alternatively,  it  may  be  chosen  to  coincide  with  the  walls  of  a  wind  tunnel  or 
towing  tank.  In  the  girthwise  direction,  NN  stations  are  used,  with  s  =  2  and 
NN  -  1  corresponding  to  the  keel  (0  =  0°)  and  the  waterplane  (0  =  90°), 
respectively;  and  c  =  1  and  NN  are  used  to  enforce  the  plane-of-symmetry 
conditions.  In  order  to  avoid  a  sudden  change  of  numerical  grid  near  the 
stern  region  where  the  hull  ends  abruptly,  an  imaginary  wake  centerplane  is 
specified  to  ensure  a  smooth  variation  of  geometric  coefficients  in  this 
region.  The  vertical  extent  of  the  wake  centerplane  is  reduced  linearly  up  to 
a  station  5  =  LE  in  the  wake.  The  depth  of  the  centerplane  is  kept  constant 
beyond  that  station. 

In  general,  the  grid  control  functions  f1  are  three-dimensional  and  can 

be  arbitrarily  specified  to  yield  the  desired  grid  distribution.  There  are, 

however,  no  general  rules  for  the  determination  of  the  most  appropriate  grid- 

control  functions.  For  the  present  study,  the  grid-control  functions  use!  in 

Chen  and  Patel  (1985b)  were  modified  and  generalized  to  provide  control  of  the 

numerical  coordinates  in  the  radial  as  well  as  the  girthwise  directions.  In 

o 

particular,  the  function  f  which  controls  the  grid  concentration  in  then- 
direction  is  continuously  adjusted  to  achieve  a  direct  control  of  the  grid 
spacing  between  the  constant-n  lines.  The  grid-generation  technique  is  out¬ 
lined  below. 
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The  choice  of  x  =  x(5)  simplifies  the  control  function  a  because  it  is 
uniquely  determined  from  equation  (35a),  i.e., 


2a 


=  7^-  =  fn(£  only) 


(37) 


In  other  words,  the  function  a  is  related  to  the  distribution  of  the  axial 
stations,  which  is  chosen  to  concentrate  grid  points  near  the  stern  and  in  the 
near  wake. 


In  a  similar  manner,  the  function  c  is  related  to  the  grid  distribution 
in  the  circumferential  or  girthwise  direction,  6  =  tan“^(y/z).  In  general,  c 
can  be  a  function  of  C,n  and  ?  to  yield  desired  grid  spacings  between  the  con¬ 
stant-?  lines  in  the  ?-  direction.  For  the  present  applications,  however,  it 
is  sufficient  to  employ  a  one-dimensional  control  function  c(?)  which  is  fixed 
in  the  axial  as  well  as  the  radial  directions,  i.e.. 


e 

2c  =  =  fn(?  only) 


(38) 


with  0  =  e{?)  prescribed  on  the  outer  boundary  at  the  upstream  station. 

The  specification  of  the  function  b,  which  controls  the  grid  distribution 
in  the  n-direction,  requires  greater  care  since  it  must  satisfy  several  con¬ 
flicting  requirements  at  the  same  time.  First,  in  the  present  treatment  of 
the  wall  boundary  conditions  using  wall  functions,  it  is  necessary  to  require 
at  least  two  near-wall  grid  points  to  lie  in  the  logarithmic  law-of-the-wal 1 
region  (say  50  <  y+  <  500).  Second,  it  is  desirable  (although  not  necessary) 
to  have  an  orthogonal  grid  in  the  wall  region  to  facilitate  the  application  of 
the  boundary  conditions.  Third,  the  grid  concentration  must  be  such  that 
there  is  a  sufficient  number  of  points  across  the  boundary  layer  whose  thick¬ 
ness  varies  greatly  in  the  axial  as  well  as  in  the  girthwise  directions. 
Fourth,  and  perhaps  the  most  difficult,  is  the  problem  of  obtaining  a  suffici¬ 
ently  accurate  solution  in  the  neighborhood  of  geometric  singularities.  It  is 
obvious  that  a  three-dimensional  control  function  is  needed  to  make  the  local 
adjustments  to  meet  all  these  requirements. 

As  noted  above,  the  grid-control  functions  used  in  Chen  and  Patel  (1985b) 
were  modified  by  adding  a  three-dimensional  local  correction  function  to 
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provide  greater  control  of  the  grid  distribution  near  the  hull  surface  and  the 
wake  centerplane.  We  choose 


b  =  b  +  b 
o  c 


(39) 


where  bQ  is  the  two-dimensional  grid  control  function  given  in  Chen  and  Patel 
(1985b),  i.e., 


2b, 


1  a33 

=  t  (rnn+  *22  hi'  22 
n  9  rg 


=  2 


i  „33 

z  znn  22  z«  ’ 
*  9 


=  fn(S,n) 


(40) 


Since  this  control  function  is  obtained  from  the  prescribed  grid  distribution 
z  =  z(5,n)  and  the  associated  transverse  curvature  z??  on  the  keel  plane  C  = 
2,  there  is  very  little  control  for  the  grid  distribution  away  from  the 
keel.  This  was  evident  from  the  results  presented  in  Chen  and  Patel  (1985b) 
which  showed  a  lack  of  near-wall  grid  concentration  around  the  waterplane, 
especially  in  the  stern  and  the  near-wake  regions.  In  order  to  improve  the 
resolution  in  the  wall  layer  and  the  near  wake,  a  three-dimensional  correction 
function  is  applied  to  control  the  normal  distance  distribution  in  the  girth- 
wise  direction  by 

bc  =  cA(?,t)e"d(n'2>  (41) 
wbere  M5  ,s)  .  , 

and  n(5,2,C)  is  the  normal  distance  between  n  =  2  and  the  body  surface  which 
is  n  =  1.  The  control  parameters  c  and  d  are  chosen  to  yield  the  desired  grid 
distributions.  Note  that  A(5,?)  =  0  if  n(5 ,2 ,c)  =  n(?,2,2)  for  all  girthwise 
stations.  In  other  words,  if  the  control  function  b  is  updated  continuously, 
the  final  correction  function  will  be  zero  when  the  solution  converges,  and 
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the  normal  distance  thus  obtained  will  be  constant  in  the  girthwise  direc¬ 
tion.  It  is,  however,  not  necessary  to  satisfy  this  relation  exactly  although 
it  is  desirable  to  be  able  to  control  the  near-wall  grid  spacing  directly. 
For  the  present  calculations,  a  compromise  solution  was  arrived  at  by  updating 
b  only  a  few  times  during  the  iterative  process  and  then  fixing  it  when  a 
satisfactory  grid  distribution  was  reached.  The  correction  function  bc  was 
updated  every  20  iterations  between  iteration  numbers  40  and  200,  and  then 
kept  unchanged  for  subsequent  iterations.  The  grid  control  function  thus 
obtained  was  enforced  only  upto  a  station  £  =  MIB  near  the  stern.  Beyond 
that,  the  f^  function  is  reduced  linearly  to  zero,  and  remained  zero  in  the 
far  wake.  For  consistency,  the  grid  distribution  z(S,n,2l  on  uie  keel  plane 
was  allowed  to  move  freely  based  on  the  prescribed  function  b  in  this 
region.  It  should  be  noted  here  also  that  the  correction  function  diminishes 
exponentially  away  from  the  hull  surface  and  wake  centerplane.  This  leads  to 
a  nearly-orthogonal  numerical  grid  for  much  of  the  solution  domain  while 
providing  the  control  of  the  grid  spacing  in  the  wall  layer. 

Equations  (35),  with  the  control  functions  specified  by  equations  (37) 
through  (41),  were  solved  by  an  exponential-linear  scheme  described  in  Chen 
and  Patel  (1985a).  With  this  scheme,  the  discretized  equations  (35)  can  be 
written 


(2g^a  coth  a  +  2g2^b  coth  b  +  2g^c  coth  cL  _<f>F 

^  s  j'l  js 

*  <gna  csch  a)S(lw5<e'\+w+  e\.w> 
t  (g22b  csch  f)Cjn>t(e'b+5>I,tlj?+ 

+  (g33c  csch  c)5ini5(e-c+?ini5+f 

+  0,5ge! n,C(<f,?+l,n+l,?+  Vl,n-l,c"  *5+l,n-l,c‘  *5-l,n+l,c) 


+  0,5g  5,n,c^5+l,n,c+l+  ^C-l.n.c-l'  ^S+l.n  ,c-l~  <*>?-l ,n  ,c+l  ^ 


+  0.5gf„  .(» 


/  m  +  <J>  «•  (i)  —  d)  ) 

5,n,5'9C,n+l,c+l  C,n-l,c-l  C,n+1,4-1  5,n-l,c+l 
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where  <f  represents  either  y  or  z.  With  the  Neumann  boundary  conditions  speci¬ 
fied  on  all  boundaries  of  the  constant-5 1  stations,  equations  (42)  are  solved 
by  a  tridiagonal-matrix  algorithm  using  line-by-line  iterations. 


(b)  Finite-Analytic  Discretization  for  the  Transport  Equations 

In  order  to  correctly  handle  the  elliptic  nature  of  the  flow  past  arbi¬ 
trary  three-dimensional  shapes,  the  finite-analytic  method  of  Chen  and  Chen 
(1982,  1984)  is  revised  and  extended  to  solve  the  five  transport  equations  for 
mean  velocities  (U,V,W)  and  turbulence  quantities  (k,e).  The  most  general 
version  of  this  method  would  involve  an  analytic  solution  of  the  linearized 
transport  equations  in  a  three-dimensional  element  and  would  result  in  a  28- 
point  discretization  formula.  While  such  a  scheme  may  be  required  for  the 
solution  of  highly  three-dimensional  flows,  for  applications  to  many  aerody¬ 
namic  and  hydrcJvnamic  problems  it  suffices  to  use  a  simplified  method  to 
reduce  computer  time  and  storage.  Here,  we  adopt  a  hybrid  method  which  com¬ 
bines  a  two-dimensional  local  analytic  solution  in  the  nc-plane  with  a  one¬ 
dimensional  local  analytic  solution  in  the  5-  direction.  Details  of  this 
numerical  scheme  are  described  in  the  following. 

In  the  finite-analytic  approach,  equations  (24)  are  locally  linearized  in 
each  rectangular  numerical  element.  A?  =  An  =  a?  =  1,  by  evaluating  the  coef¬ 
ficients  of  the  convective  terms  at  the  interior  node  P  of  each  local  element 
(Figure  4),  i.e., 

9p'*CC+  V*nn+  2<VpV  2<VpV  2(VpV  <VpV  (Vp  (43) 


Introduction  of  the  coordinate-stretching  functions 


* 

5 


* 


n 


n 


(44) 


in  equation  (43)  reduces  it  to  the  standard  three-dimensional  convective- 
transport  equation  described  in  Chen  and  Chen  (1982,  1984),  i.e., 


+  A  +  A 
Y  *  Y  *  * 


=  2C4>  *+  2B4>  *+  2A<j>  *+ 


5  K 


n  n 


4  ? 


RV 


<Vp 


(45) 
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where  A  = 


(AJ 


li£ 


,  33 
/gP 


B  = 


(Vp 


22 


c  = 


((Vp 

,“11 

/gP 


>  R  =  <Vp 


for  a  numerical  element  with  dimensions 


*  1 
ac  =  i  =  — - 


•'g 


An  =  k  = 


,12 


*  1 
AC  =  h  = 

,13 

% 


(46) 


For  the  type  of  applications  considered  here,  it  is  convenient  to  decom¬ 
pose  equation  (45)  into  a  one-dimensional  and  a  two-dimensional  partial  dif¬ 
ferential  equation  in  the  following  way: 


2W  *-<(>**+  R^t  +  2<j)=  G(£  >  n  *  ^  »  t) 

?  £  £ 

•k  -k  -k 

*  *+  4>  *  *-  2B4>  *-  2A4>  *  =  G(5  ,  n  ,  C  ,  t) 
n  n  c  C  n  4 


(47) 

(48) 


If  we  further  require  the  source  functions  G  and  S,  to  be  constant  in  each 

<p 

local  element  and  the  time  derivatives  to  be  approximated  by  a  backward-dif¬ 
ference  formula,  equations  (47)  and  (48)  reduce  to  the  standard  one-  and  two- 
dimensional  convective-transport  equations  described  in  Chen  and  Chen  (1982), 
respectively.  The  analytic  solution  of  the  one-dimensional  equation  (47)  can 
be  readily  obtained  as: 

★ 

4>  =  a(e2C?  -1)  +  bC*+  c  (49) 


By  substituting  the  exponential-linear  solution  (49)  into  equation  (47),  the 
source  function  G(0, 0,0,0)  =  g  becomes 


where  the  subscripts  U  and  0  denotes  the  upstream  and  downstream  nodal  values 
(Figure  4),  respectively,  superscript  (n-1)  denotes  the  value  at  the  previous 
time  step,  and  t  is  the  time  step. 

By  specifying  a  combination  of  exponential  and  linear  boundary  functions, 
which  are  derived  from  the  natural  solutions  of  the  governing  equations,  on 
all  four  boundaries,  o  =  ±  k  and  ?  =  +  h,  of  the  transverse  section  of  each 
local  element  (5-plane),  i.e., 


♦(k,c*)  =  an(e2AC  -  1)  +  bnt*+  cn 
=  a s(e2Ac  -  1)  +  b$c*+  c$ 

*  (51) 

4>{n  ,h)  =  ae(e  +  bgn  +  cg 

'k 

<f> (n*,-h)  =  aw(e2Bn  -  1)  +  bwn*+  cw 

where  a,  b  and  c  are  constants,  the  two-dimensional  equation  (48)  can  be 
solved  analytically  by  the  method  of  separation  of  variables  or  any  other 
analytic  technique.  Details  of  the  solution  procedure  are  described  in  Chen  & 
Chen  (1982,  1984).  When  the  local  analytic  solution  thus  derived  is  evaluated 
at  the  central  node  P  of  the  element,  the  following  nine-point  finite-analytic 
algebraic  equation  is  obtained: 

^P  =  CNE*NE+  WnW+  CSE*SE+  Sw^SW  +  CEC4>EC+  VwC 

+  Vnc*  CSC*SC~  cp9  (52) 


30 


where 


CSC  =  ^2cosh  Bk5  PA 


C  =  e"2BkC 
UNC  e  LSC 


r  =  f - § - )  p 

cWC  11 2  cosh  Ah'  B 

C  =  e~2^C 
lEC  e  tWC 

e  Ah+Bk 

CSW  =  *-4  cosh  Ah  cosh  Bk-*  (1“  PA“  PB) 


r  —  p“2Ahp 
CSE  "  e  CSW 


CM,,  =  e‘"KCc 


r  _  p-2Ah-2Bk r 
CNE  "  e 


r  h  tanh  Ah  ,,  D  .  k  tanh  Bk  n  D  « 
Cp  =  2A  U-PA)  -  2B  U-HgJ 


P A  =  4E2  Ah  cosh  Ah  cosh  Bk  coth  Ah 

p  _  i  j.  Bh_coth_Bk  .  _ 

^B  Ak  coth  Ah  vlA  u 


Co  *  i’ 


-  (-l)m(Xmh) 


™'1  [(Ah)2+  (Xmh)2]2cosh  /A2+  B2+  >2  k 


X  n  =  (in  -  -k)  * 
m  2 

The  above  coefficients  are  simple  rearrangements  of  those  given  in  Chen  & 
Che  „  but  are  mo'e  convenient  for  efficient  numerical  calculations.  Note  that 
icme  of  tnem  are  interrelated  as  follows: 


CNC  +  CSC  "  PA 


CEC  +  CWC  =  PB 


(53) 


CNE  +  CNW  +  CSE  +  CSW  =  1  ‘  PA  '  PB 

Since  both  PA  and  Pg  are  positive  and  (Pft  +  Pg)  <  1,  the  finite-analytic 
coefficients  are  always  positive.  Also,  the  exponents  in  these  coefficients 
provide  a  gradual  upwind  bias  as  the  cell  Reynolds  numbers  |2Ah|  and/or  |2Bk| 
increase.  Thus,  the  behavior  of  the  convective-diffusion  equation  is  properly 
captured  and  numerical  diffusion  is  minimized  due  to  the  inclusion  of  all 
corner  points.  For  large  cell  Reynolds  numbers,  the  series  summation  in  E£ 
can  $e  avoided  by  considering  the  asymptotic  expressions  of  PA  and  Pg  based  on 
the  theory  of  characteristics,  i.e., 

Ak  coth  Ah  >  Bh  coth  Bk  :  Pft=  0,  PD=  1  -  Bh  coth  Bk/Ak  coth  Ah 

A  d  (54) 

Ak  coth  Ah  <  Bh  coth  Bk  :  Pg=  0,  P^=  1  -  Ak  coth  Ah/Bh  coth  Bk 

Since  the  "downstream"  influence  is  negligible  at  large  cell  Reynolds  numbers, 
the  above  approximations  do  not  introduce  a  significant  error  in  the  solution, 
but  the  computing  time  is  greatly  reduced. 

By  substituting  the  nonhomogeneous  term  g  from  equation  (50)  into  equa¬ 
tion  (52),  a  twelve-point  finite-analytic  formula  for  unsteady,  three-dimen¬ 
sional,  elliptic  equations  can  be  obtained  in  the  form 


*P  i  4.  r  rr  4.  r  4.  &  1  {CNE*NE+  CNVJ*NW+  CSE*SE+  CSW<f>SW+  CEC4,EC+  VwC 
I  +  tpLLu+  Ld+  t  J 


R  ,  n-1 


+  CNC<f>NC+  CSC^SC+  Wu+  CD*D+  t  ^P  *  ”  Cp(Vp} 


(55) 


or 


8 


3 - 5  C.  ♦  .  +  CD(C,rfll+  CJn+  J^"1) 


l+Cp[Cu+CD+?]  1  nb  nb  P'  U  U  D  D  T  P 


(55a) 
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where  the  subscript  nb  denotes  neighboring  nodes  (NE:northeast,  NW:northwest, 

etc.).  It  is  seen  that  <j>p  depends  upon  all  the  eight  neighboring  nodal  values 

in  the  transverse  plane  as  well  as  the  values  at  the  upstream  and  downstream 

nodes  ^  and  and  the  values  at  the  previous  time  step  <f>p  .  When  the  cell 

Reynolds  number  2C  becomes  large,  C..  2C/£  =  (C.  )n  and  C.,  0,  and  equation 

U  r  U 

(5 5)  reduces  to  the  partially-parabolic  formulation  of  Chen  and  Patel  (1985a) 

which  used  =  (C^)p  and  CQ  =  0.  Thus,  the  extension  of  the  earlier  par¬ 

tially-parabolic  method  to  the  present  ful ly-elliptic  form  is  straightforward, 
and  if  <J>D  is  obtained  from  the  previous  time  step,  the  same  algorithm  can  be 
employed  for  both  formulations. 

Since  equations  (55)  are  implicit,  both  in  space  and  time,  at  the  current 

station  of  calculation,  their  assembly  for  all  elements  results  in  a  set  of 

simultaneous  algebraic  equations.  These  equations  are  solved  by  the  tridiag¬ 
onal-matrix  algorithm.  Because  it  is  not  necessary  to  obtain  a  fully- 
converged  intermediate  solution  for  steady  flows,  only  ten  line-by-line  inter¬ 
nal  iterations  are  used  during  each  global  sweep.  Furthermore,  the  finite- 
analytic  coefficients  appearing  in  equations  (55)  are  not  updated  during  these 
internal  iterations  for  economy  of  computation  time. 

(c)  Solution  of  the  Continuity  Equation:  Velocity-Pressure  Coupling 

If  the  pressure  is  known,  equations  (55)  can  be  employed  to  solve  the 
five  convective-diffusion  equations  (24)  for  U,V,W,k  and  e.  However,  the 
pressure  is  not  known  a  priori  and  must,  be  determined  by  requiring  the  velo¬ 
city  field  to  satisfy  the  equation  of  continuity  (32).  Since  a  direct  method 
for  the  simultaneous  solution  of  all  six  equations  is  not  feasible  with  pres¬ 
ent  computer  capacity,  it  is  necessary  to  convert  the  equation  of  continuity 
into  an  algorithm  for  the  calculation  of  the  pressure  field.  The  SIMPLER 

algorithm  of  Patankar  (1980)  has  been  modified  and  extended  for  this  purpose. 
A  staggered-grid  system  is  adopted.  Figure  5  shows  the  locations  of  the  nodes 
for  J,V,W,  and  p  in  this  grid.  The  turbulence  quantities  k  and  e  are  evalu¬ 
ate  at  the  pressure  nodes.  The  dashed  lines  represent  the  control  volume 
•aces,  and  the  pressure  is  calculated  at  the  center  of  the  control  volume. 
For  convenience,  'J^*  Vn,  Wg  and  pp  in  Figure  5  are  assigned  the  same  index, 
i.e.,  they  are  denoted  by  V?n?,  and  p?  ,  respectively.  Recall 
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that  the  velocity  components  U,V,  and  W  are  in  the  longitudinal,  radial  and 
circumferential  directions.  In  other  words,  they  are,  in  general,  neither 
perpendicular  to  the  control  surfaces  nor  in  the  direction  of  the  coordinate 
lines.  However,  due  to  the  choice  of  the  cylindrical  polar  coordinates  in  the 
physical  plane,  these  components  do  not  become  parallel  to  the  control  sur¬ 
faces. 

In  the  staggered  grid,  the  twelve-point  finite-analytic  formulae  (equa¬ 
tion  55)  for  the  momentum  equations  yield  the  velocity  components: 


U,  = 


V_  = 


R.  cnbunb+  +  *Vdud+  t  ud  *Vd® 

Wd(VCD+A 


1 


—  r7~«  CnbV  WA  *  <Vn» 

,+W  V  t  »„ 


1 


- R—  «  cnbwnb+  *  (CD)eWd+  ^  <ye]) 

'WV.  V 

(56) 


where  Cj,  Cn  and  Ce  are  the  finite-analytic  coefficients  Cp  evaluated  at  the 
staggered  velocity  nodes  d,  n  and  e  in  Figure  5.  Note  that  the  above  equa¬ 
tions  contain  the  pressure-gradient  terms  inside  the  source  functions.  An 
equation  for  this  unknown  pressure  field  is  obtained  as  follows. 

If  we  decompose  the  actual  velocity  field  (U,V,W)  in  the  momentum  equa¬ 
tions  (56)  into  a  pseudovelocity  field  (U,V,W)  plus  the  pressure-gradient 
terms  contained  in  the  source  functions,  i.e., 


ud =  V 


R„  {  J  (blpC+  blpn+  blp?)}d 


V„  =  V  - 
n  n 


l+Cd(Cu+  CD+  x^d 


123 
(b2Pc+  b2pn+  b2p?))n 


1+Cn(CU+  CD+  x^n 
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(57) 


C  R 

»e  *  “e - 8 - R—  {-J  (b3P5+  »k»e 

l+Ce(Cu+  c0+  T  )e 

so  that  the  pseudovelocities  contain  no  pressure  terms,  then  an  equation  for 
pressure  can  be  derived  by  requiring  the  velocity  field  to  satisfy  the  discre¬ 
tized  equation  of  continuity  (32),  i.e.. 


(dJu  +  b^V  +  b^W)^-  (bju  +  b^V  +  b3W)u+  (b'jU  +  b^V  +  bgU^ 

-  (b?U  +  blv  +  b?W)  +  (b?U  +  b^V  +  b?W )  -  (b^U  +  b?V  +  b^W)  =  0  (58) 

The  resulting  pressure  equation  will  contain  many  pressure  nodes  (see  Muraoka 
(1980,  1982)  for  example)  if  nonorthogonal  coordinates  are  employed.  It  ^s 
therefore  desirabje^  to  simplify  the  algorithm  by  introducing  modified 
pseudovelocities  (U,V,W)  by  decomposing  the  velocity  components  as  follows: 


Ud  =  Ud  '  dd(pD"  PP^ 
dn(pNC'  PP} 
We  =  We  "  de(pEC"  PP} 


where 


dd  = 


{RU  bl)dCd 


odn  ♦  cd(cu+  Cp+  ^)d] 


<RV  Wn 


•V1*  W  Vt  )„] 


(59) 


(59a) 
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The  modified  pseudovelocities  (U,V,W)  defined  above  still  contain  part  of  the 
pressure-gradient  terms  if  the  coordinate  system  is  nonorthogonal  (i.e.,  b^  * 
for  i  *  j,  see  equation  57).  These  pressure-gradient  terms  can  be  evaluated 
from  the  pressure  field  known  at  the  previous  time  step  or  iteration  without 
losing  any  accuracy  or  generality.  If  we  require  the  velocity  field  to 
satisfy  the  equation  of  continuity  (58),  a  simpler  pressure  equation  can  be 
derived  in  terms  of  the  modified  pseudovelocities  (U,V,W).  Note  that  eighteen 
velocity  components  are  involved  in  equation  (58)  for  each  control  volume. 
However,  due  to  the  staggered  grid  system  employed  here,  only  six  of  these, 
namely,  Ud,  Uu,  Vn,  V$,  W  ,  and  W  ,  can  be  obtained  directly  from  the  govern¬ 
ing  equations  (56).  It  is,  therefore,  necessary  to  approximate  the  remaining 
twelve  by  interpolations.  A  simple  linear  interpolation  is  used  here  to 
evaluate  these  from  the  velocity  field  known  at  the  previous  time  step  or 
iteration,  so  that  the  continuity  equation  becomes 

(b]u)d  -  (b]u)u+  (b^ V) n~  <b2V)s+  (b3W)e'  (b3W)w+  D1  =  0  {60) 

where 

D1  =  (b2V  +  b3W)d  ‘  (b2V  +  b3W)u+  (blU  +  b3W)n'  (blU  +  b3W)s 

+  (b^U  +  b^V)  -  (b^U  +  b~V)  (60a) 

I  2  e  I  2  w 

is  the  mass  source  obtained  from  the  velocity  field  at  the  previous  time  step 
or  sweep.  An  equation  for  pressure  is  then  derived  by  substituting  equations 
(59)  into  (60),  i.e., 

apPp  =  adpD  +  auPy  +  anpNC  +  aspsc  +  aepEC  +  awpwc  -  D  (61) 


where 


au  ■  <bl>u  du 
an  "  dn 

as  =  (b2}s  ds  (61a) 

ae  '  <b3>e  de 
aw  '  (b3*w  dw 

aP  =  ad  +  au  +  ae  +  aw  +  an  +  as 
and 

D  =  (b]  U)d-  (b]  U)u+  (b|  V)n-  (b^  V)s  +  (b|  W)e-  (b|  W)w-  D1  (61b) 

AAA 

The  modified  pseudovelocities  (U,  V,  W)  contain  the  neighboring  nodal 
values  of  velocity,  source  functions,  and  part  of  the  pressure-gradient  terms. 
All  of  them  can  be  evaluated  from  the  information  known  at  the  previous  time 
step  or  iteration.  Therefore,  apart  from  the  interpolations  for  D] ,  the  pres¬ 
sure  equation  (61)  is  still  an  exact  algebraic  representation  of  the  equation 
of  continuity  (32).  In  this  fashion,  the  pressure  field  can  be  obtained 
directly  from  an  estimated  velocity  field. 

Although  the  guessed  pressure  field  can  be  updated  directly  by  equation 
(61),  in  practice  the  new  pressure  field  may  produce  a  velocity  field  which 
does  not  satisfy  the  equation  of  continuity.  An  iterative  procedure  is  there¬ 
fore  required  to  correct  this  erroneous  velocity  field  to  achieve  more  rapid 
convergence.  Here,  a  velocity-correction  formula,  similar  to  that  used  in  the 
SIMPLE  algorithm,  is  derived  in  terms  of  the  pressure-corrections. 

If  we  denote  the  imperfect  velocity  field  resulting  from  an  imperfect 
pressure  field  p*  by  (U*,V*,W*),  then  the  discretized  momentum  equations  (59) 
can  be  written  as 

+  A  *  *  -k 

Ud  =  Ud  de(pD"  M 
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(62) 


l 


k  *  ic  -k  k 

\  -  \  -  dn‘PNC-  Pp  > 

"e  =  “e*-  ‘MPEC-  Pp> 

In  order  to  improve  the  guessed  pressure  field,  such  that  (U*,V*,W*)  will 
eventually  satisfy  the  equation  of  continuity,  one  needs  to  know  how  the 
velocity  components  respond  to  a  change  in  the  pressure  field.  Such  a  rela¬ 
tion  can  be  obtained  by  subtracting  equation  (62)  from  equation  (59),  i.e., 

ud  -  ud  =  <V  •  dd<PD  -  Pp> 

V„  -  \  ■  <V  V  -  dn(PNC-  Pp>  (63) 

“e  -  “e  "  <"e-  “e*>  *  ^'4-  4 

where  p'  =  p  -  p*  is  the  pressure  correction,  and  (U-U*),  (V-V*)  and  (W-W*) 

are  the  corresponding  velocity  corrections.  If  we  require  the  velocity  field 
to  satisfy  the  equation  of  continuity  (60),  an  equation  for  the  pressure 
correction  p'  can  be  derived.  However,  due  to  the  implicit  nature  of  the 
velocity  corrections  arising  from  the  pseudovelocities,  the  resulting  pressure 
correction  equation  would  involve  the  pressure  corrections  at  all  grid  points. 
It  is  not  necessary  to  retain  such  a  complicated  formulation  because  both  the 
pressure-  and  velocity-corrections  are  zero  in  the  final  converged  solution. 
Since  both  the  pressure-  and  velocity-corrections  become  trivial  when  the 
solution  converges,  it  is  possible  to  omit  that  part  of  the  velocity-correc- 

A  A  ^  A  A  ^  A  A  ^ 

tions,  (U  -  U  ),  (V-V  )  and  (W  -  W  ),  which  represents  the  indirect  influ¬ 
ence  of  velocity  corrections.  With  this  approximation,  the  velocity-correc¬ 
tions  are  expressed  explicitly  in  terms  of  the  pressure-corrections  as 


Ud  =  Ud  "  dd*pD“  PP* 
dn{pNC-  PP> 
We  =  We  “  de(pEC‘  PP* 
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By  requiring  the  velocity  field  to  satisfy  the  equation  of  continuity  (58),  a 
pressure-correction  equation  is  then  obtained  in  the  form 

aPpP=  adpD+  aupU+  anpNC+  aspSC+  aepEC+  awpWC~  D  (65) 

with 

D*=  (b]u*)d  -  (b]u*)u+  (b^V* )n~  (b^V* ) s+  (b|w*)e~  (b|w* )w«  D]  (65a) 

where  the  coefficients  ap,  ad,  etc.,  are  as  defined  in  equations  (61a).  Note 
that  the  pressure-correction  equation  (65)  is  similar  to  the  pressure  equation 
(61).  Although,  unlike  the  pressure  equation,  the  pressure-correction  equa¬ 
tion  is  not  exact,  the  approximations  made  influence  only  the  rate  of  conver¬ 
gence  but  not  the  final  converged  solution. 

The  systems  of  algebraic  equations  formed  by  the  assembly  of  the  pressure 
and  pressure-correction  equations,  (61)  and  (65),  respectively,  are  solved  by 
the  tridiagonal-matrix  algorithm  with  several  line-by-line  internal  iter¬ 
ations.  The  finite-analytic  coefficients  ap,  ad,  etc.,  are  updated  in  each 
upstream  to  downstream  global  sweep,  but  remain  the  same  during  the  internal 
iterations. 

(d)  Boundary  Conditions 

For  the  calculation  of  ship  stern  and  wake  flows  considered  in  this 
report,  it  is  assumed  that  the  ship  is  symmetric  about  a  vertical  centerplane 
and  the  waterplane  is  regarded  as  a  plane  of  symmetry.  The  appropriate  bound¬ 
ary  conditions  are  then  as  follows  (see  Figure  2). 

(1)  Initial  or  Upstream  Section  (C  =  1) 

The  distribution  of  the  velocity  components  (U,V,W)  and  the  turbulence 
parameters  (k,e)  are  assumed  known  at  an  upstream  transverse  section  either 
from  detailed  boundary-layer  calculations  or  from  simple  correlations.  Bound¬ 
ary  conditions  for  pressure  and  pressure-correction  are  not  required  in  the 
present  staggered-grid  arrangement  since  the  pressure  is  implicitly  determined 
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by  the  specification  of  the  velocity  components.  Also,  since  the  velocity 
components  are  known,  there  are  no  corrections  to  be  made,  i.e.  (U,V,W)U  = 
(U*,v\w*)u  =  (U,V,W)U. 

(2)  Exit  or  Downstream  Plane 

The  exit  plane  is  usually  located  in  the  wake  far  downstream  from  the 
stern,  and  the  zero  pressure  gradient  (p^=  0)  condition  is  specified  there. 
Since  axial  diffusion  in  the  far  wake  is  negligible,  the  exit  conditions  for 
the  transport  quantities  <t>  =  U,V,W,k,e  in  the  momentum  and  turbulence-model 
equations  are  simply  <f>^  =  0. 

(3)  Body  Surface  (n=l) 

For  laminar  flows,  the  numerical  solution  is  usually  carried  out  upto  the 
solid  surface  using  no-slip  conditions,  i.e.,  U=V=W=Q  (see,  for  example,  Chen 
and  Patel,  1987a).  Strictly  speaking,  the  same  conditions  should  also  be  used 
in  turbulent  flow  calculations.  Tnis,  however,  will  require  a  large  number  of 
grid  points  to  resolve  the  large  gradients  in  the  near-wall  region  and,  more 
importantly,  an  appropriate  near-wall  turbulence  model  to  account  for  the 
wall-proximity  effects.  In  view  of  the  complexity  involved  in  resolving  the 
near-wall  flow,  it  is  preferable  to  employ  a  simpler  wall-function  approach 
which  avoids  the  solution  of  the  equations  of  motion  and  turbulence  model  in 
the  wall  region.  In  the  present  study,  the  two-point  wall-function  approach 
of  Chen  and  Patel  ( 1985a, b,  1988)  is  employed  with  necessary  modifications  to 
determine  the  boundary  conditions  on  a  fictitious  boundary  (n  =2 )  located  in 
the  ful ly-turbulent  logarithmic  layer. 

The  present  wall-function  approach  differs  from  the  usual  practice  in  two 
respects.  First,  the  effects  of  pressure  gradients  on  the  flow  in  the  wall 
region  are  taken  into  account  by  the  use  of  a  generalized  law  of  the  wall 
given  by  Patel  (1973),  i.e., 

„  1  r4  n  +\/)1/2  -  i  +  1/2 

S-  =  7  nntj - — TTp - ]+2[(l  +  A TyV  -  1]}  +  B  +  3.7  A  (C6) 

UT  \  (1  +  A  v+)!/2+  1  T  p 
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¥ 


2  + 

in  which  U  is  the  friction  velocity  defined  by  U  =  /t  /pu  ,  y  =  Re  U  y 

T  TWO  T  J 

is  the  dimensionless  distance  measured  in  the  direction  normal  to  the  surface, 
3 

Ap  =  Vp/Re(JT  is  the  dimensionless  pressure  gradient,  at  is  the  dimensionless 
shear-stress  gradient  which  is  approximated  by  1/2  A  ,  q  is  the  magnitude  of 
the  velocity,  <  =  0.42  is  the  von  Karman  constant,  and  B  =  5.45.  Second,  we 
ensure  that  at  least  two  points  (n  =  2  and  3  in  Figure  6)  are  located  in  the 
logarithmic  region  and  explicitly  satisfy  equation  (66)  at  both.  This  avoids 
the  need  for  a  separate  analysis  for  the  flow  between  the  wall  and  the  first 
near-wall  mesh  point  which  is  used  in  almost  all  previous  applications  of  wall 
functions.  For  the  present  calculations,  as  in  the  earlier  study  of  Chen  and 
Patel  (1985a, b,  1988),  equation  (66)  is  employed  to  improve  the  prediction  of 
the  wall  shear  stress  and  the  associated  boundary  conditions  in  adverse  pres¬ 
sure  gradients,  while  the  usual  logarithmic  law  of  the  wall  (i.e.,  a^  =  Ap  = 
0)  is  used  in  favorable  pressure  gradients. 

In  the  present  procedure,  a  value  for  U  is  assumed  and  the  boundary 
conditions  at  n  =  2  are  determined  from  equation  (66)  and  the  assumptions  of 
local  equilibrium  for  k  and  e  ,  i.e.. 


1/2 


9?  ]  4  (1  +  ^-^9)  “  f  +  1/9 

r  -  <  7: — I-h7rT]+2[<i +  V2> 7  - 1]> +  B  +  3-7  A, 

t  x  (1  +  ATy2)  +  1  ' 


4 

2  /c 


(67a) 


(67b) 


e2  =  <y. 


(67c) 


The  numerical  solution  then  provides  the  velocity  at  n  =  3  and  is  updated 
by  requiring  this  velocity  also  to  satisfy  equation  (66),  i.e.,  by  solving 


qo  -1  4  d+W  -1  +17? 

TT  =  ^nn[J - - ]+2C(l  +  i.y  )  '  -  1]>  +  B  +  3.7  A  (68) 

UT  T  (1  +  V3}  +  1  P 
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by  a  root-finding  technique.  Thus,  an  iterative  procedure  is  used  to  satisfy 
the  wail  boundary  conditions  in  the  case  of  turbulent  flow,  and  to  deter¬ 
mine  UT.  Typically,  five  iterations  are  required  to  obtain  satisfactory 
convergence. 

By  anchoring  the  solution  at  two  near-wall  mesh  points  on  the  law  of  the 
wall,  the  present  wall-function  approach  removes  much  of  the  sensitivity  of 
the  numerical  solution  to  the  location  of  the  first  mesh  point  which  has  been 
observed  in  previous  treatments.  The  procedure  Is  quite  straightforward  for 
two-dimensional  and  axisynmetric  flows.  For  three-dimensional  flows,  however, 
an  additional  assumption  concerning  the  direction  of  the  velocity  vector  is 
required  to  determine  the  individual  components  (U,V,W)  because  the  law  of  the 
wall  gives  only  the  velocity  magnitude.  In  Chen  and  Patel  (1985a),  it  was 
assumed  that  q2  is  parallel  to  the  wall  and  there  is  no  rotation  of  the  veloc¬ 
ity  vector  between  n=2  and  n=3  in  planes  parallel  to  the  surface.  Although 
the  latter  approximation  resulted  in  some  simplification  in  the  application  of 
the  wall-function  boundary  conditions,  it  led  to  the  underprediction  of  the 
secondary  velocities  in  the  stern  and  wake  regions.  The  approach  of  Chen  and 
Patel  (1985a)  was,  therefore,  modified  to  account  for  the  rotation  of  the 
velocity  vector  in  the  wall  layer. 

To  use  the  law-of-the-wall  formula  (66),  it  is  convenient  to  relate  the 

velocity  components  (U,V,W)  in  the  governing  equations  to  the  physical  compon- 
C  n  ? 

ents  (q  ,  q  ,  q  )  along  the  body-fitted  coordinates  as  follows: 

a  =  (U,V,W)  =  q  t?+  qn  1^+  q  t?  (69) 

where  t  ,  t  and  t  are  unit  tangents  in  the  body-fitted  coordinates  given  by 

I5  =  (x5,  r?,  r8?) 

/g11 

In  =7=<V  V  rV  l70) 

'322 
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Thus, 


U  = 


V  = 


W  = 


re 


’ll 


and  the  inverse  relations  are 
c  *^9 1 1  I  1  1 

q^  =  -j11  (b-jU  +  b^V  +  b-Jw) 
qn  =  (b'jU  +  b^V  +  bgW) 


(b^U  +  b^V  +  b^W) 


(71) 


(72) 


The  geometric  coefficients  appearing  in  the  above  relations  are  evaluated  at 

the  wall,  n=l,  so  that  the  rotation  of  the  velocity  vector  can  be  determined 

relative  to  the  surface  coordinates.  The  projection  of  the  velocity  vector  on 
*  n 

the  surface,  q  (i.e.,  without  q  component)  and  the  angle  B  between  this  and 
the  S-direction  of  the  surface  coordinates  are  given  by 

q*  =  q%  +  q?ir  (73) 

/v  J 

tan  B  =  -fll-ilO — 

q^  +  cos  Y 
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where  y  is  the  angle  between  t  and  t  .  Alternatively,  one  can  also  evaluate 

£  r  * 

the  components  q  and  q  when  the  magnitude  and  direction  of  the  vector  q  are 
specified,  i.e., 


* 

C  _  Q  sin  6 
q  ~  sin  y 


q  =  q  cos  3  -  q  cos  Y 


(74) 


In  the  present  procedure,  we  first  impose  the  constraint  that  q2  is 
parallel  to  the  wall  and  therefore  require  q^=  0.  However,  instead  of  the  no 
rotation  (i.e.,  =  8^)  assumption  imposed  in  Chen  and  Patel  (1985a),  the 
effect  of  the  velocity  vector  rotation  in  the  near  wall  region  is  taken  into 
account  by  determining  the  direction  of  ^  using  a  quadratic  extrapolation  of 
the  angles  8^,  3^  and  3,.  of  the  velocity  vectors  at  n  =  3,  4  and  5.  The 
angle  6^  thus  obtained  is  then  used  in  equation  (74)  to  calculate  the  velocity 
components  at  n  =  2  from: 


q2  =  q2  cos  -  q^  cos  Y 

q2  Sin  8? 
q2  sin  y 


(75) 


92  =  0 

Note  that  q2  =  q^  since  it  is  assumed  that  q2  is  parallel  to  the  wall,  i.e.  q1^ 
=  0. 

The  iterative,  two-point,  wall-function  approach  for  a  three-dimensional 
flow  can  be  summarized  as  follows.  With  an  assumed  U  ,  and  therefore  q2  from 
equation  (67),  equations  (75)  give  the  components  q|,  q?,,  q^,  and  equations 
(71)  give  the  components  (l^jV^W.?)  which  are  required  as  the  boundary  condi- 


hence  93.  A  new  value  of  is  obtained  from  equation  (68),  and  equations 
(72)  and  (76)  are  used  to  find  q  and  3  at  n  =  3,  4  and  5  for  the  quadratic 
extrapolation  to  determine  8  .  Finally,  q2  is  obtained  from  equation  (67) 
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using  the  new  UT ,  q^, 
updated  boundary  conditions  for  the  numerical  solution  and  the  procedure  is 
iterated  to  convergence. 


q^,  q!)  found  from  equations  (75).  These  provide  the 


It  should  also  be  noted  that  in  evaluating  the  distance  y  normal  to  the 
wall,  required  in  the  law  of  the  wall,  it  is  necessary  to  determine  the  direc¬ 
tion  cosine  between  £n  and  jQn.  This  is  given  by 


cos  a 


2  2? 

b-,x„  +  b0r  +  b,r9„ 
In  2  n  3  n 


'99229 
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(76) 


in  which  all  the  geometric  coefficients  are  again  evaluated  at  the  wail,  i.e., 
at  n  =  1. 


(4)  Outer  Boundary 

In  order  to  simulate  the  flow  past  a  body  in  an  unrestricted  uniform 
stream,  the  outer  boundary  of  the  solution  domain  is  placed  at  a  large  dis¬ 
tance  from  the  body.  The  boundary  conditions  then  become 

U  =  l,  W  =  0,  |Hf=0,  p  =  0  (77) 

The  radial  component  of  velocity  (V)  at  the  outer  boundary  is  not  specified 
but  is  determined  by  the  solution.  The  outer  boundary  can  also  be  made  to 
coincide  with  the  walls  of  a  wind  tunnel  or  towing  tank,  in  which  case  the 
boundary  conditions  are  identical  with  those  described  in  the  previous  sec¬ 
tion. 

While  the  calculations  presented  in  this  report  are  of  the  type  described 
above,  we  also  recognize  that  the  outer  boundary  can  be  made  to  coincide  with 
a  surface  in  the  inviscid  flow  provided  appropriate  match  conditions  are 
provided.  This  viscous-inviscid  matching  approach  is  not  pursued  further. 

(5)  Symmetry  Planes 

On  the  centerplane,  i.e.,  the  ship  keel  and  the  wake  centerplane. 
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and  on  the  undisturbed  free  surface,  i.e.,  the  waterpiane,  symmetry  conditions 


W  =  0, 


DJU  2!  _  3Jk  _  9e_  _ 
3?  =  3?  =  3?  =  3? 


are  enforced  for  the  transport  quantities. 


(78) 


(6)  Initial  Conditions  (t  =  0) 

For  the  steady-flow  calculations  considered  here,  the  initial  conditions 
for  the  transport  quantities  4>  =  'l^VjWjk.e)  are  taken  directly  from  the 
values  known  at  the  immediately  upstream  station  during  the  first  sweep  and 
the  pressure  is  assumed  to  be  zero  throughout  the  flow  field.  Although  it  is 
possible  to  assume  more  realistic  initial  conditions  and  thus  accelerate  the 
convergence  of  the  solution,  these  rather  crude  initial  conditions  have  been 
used  not  only  to  simplify  the  use  of  the  method  but  also  to  demonstrate  its 
versatility. 


(e)  The  Overall  Solution  Algorithm 

With  the  grid  distributions  employed  for  the  large  Reynolds  numbers 
considered  here,  the  cell  Reynolds  number  in  the  £  -  direction  is  quite  large. 
Consequently,  in  equations  (55),  the  influence  coefficient  Cy  is  much  smaller 
than  Cy.  In  other  words,  the  streamwise  diffusion  is  small  compared  to  stream- 
wise  convection  and  transverse  diffusion.  This  enables  us  to  adopt  a  partial- 
ly-parabolic  solution  algorithm  with  minor  modifications,  rather  than  a  fully 
iterative  scheme,  to  solve  the  elliptic  equations.  As  noted  earlier,  the  pre¬ 
sent  twelve-point  finite-analytic  formula  contains  only  one  additional  influ¬ 
ence  coefficient  Cg  which  is  not  present  in  a  partially-parabolic  formula¬ 
tions.  If  <j>D  is  evaluated  using  the  value  known  from  previous  sweep  or  time 
step,  then  the  partially-parabolic  solution  algorithm  of  Chen  and  Patel 
(1985a)  can  be  employed  also  to  solve  the  elliptic  equations.  Thus,  the 
influence  of  streamwise  diffusion  in  the  ful ly-elliptic  equations  is  conven¬ 
iently  taken  into  account  without  significantly  increasing  computer  storage  or 
computation  time. 

For  transient  problems,  where  the  initial  and  boundary  conditions  are 
properly  specified,  the  overall  numerical  solution  procedure  may  be  summarized 
as  follows: 
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1.  Construct  the  body-fitted  coordinate  system  for  the  given  body  shape 
and  solution  domain,  and  calculate  the  geometric  coefficients 

b'-,  g1J,  J,  etc.  from  equations  (27)  to  (31). 

2.  Specify  the  initial  conditions  for  the  velocity  and  turbulence  fields. 
Set  p  =  0  everywhere  initially. 

3.  Specify  the  velocity  and  turbulence  profiles  at  the  first,  station  5  =  1 
(these  may  be  time  dependent). 

4.  Calculate  the  finite-analytic  coefficients  for  momentum,  pressure,  and 
pressure-correction  equations  at  the  downstream  station  from  equations 
(52)  and  (61),  respectively.  Store  only  the  finite-analytic  coeffici¬ 
ents  a^,  cr.  and  a.?  fo~  the  pressure  equation. 

5.  Solve  the  momentum  equations  ^ased  on  the  updated-pressure  field  to 
obtain  the  starred  velocity  fie  id  ( U*,  V*,  W*).  This  system  of  alge¬ 
braic  equations  is  solved  by  a  tridiagonal  matrix  algorithm. 

6.  Calculate  the  mass  source  0*,  and  solve  the  pressure  correction  equa¬ 
tion  (65)  by  tridiagonal  matrix  algorithm. 

7.  Correct  the  velocity  field  using  the  velocity-correction  formulae  (64), 
but  do  not  correct  the  pressure  field. 

8.  Update  wall-function  boundary  conditions  using  the  newly-obtained 
velocity  field  and  repeat  steps  (5)  to  (7)  for  several  internal  itera¬ 
tions. 

AAA 

9.  Calculate  the  pseudovelocities  :'U,  V,  W)  in  terms  of  the  velocity  field 
from  equation  (59).  Store  only  D  for  later  use. 

10.  Solve  equations  (55)  for  turbulence  quantities  (4>  =  k,e)  by  tridi¬ 
agonal-matrix  algorithm. 

11.  March  to  the  next  downstream  station  and  repeat  steps  (3)  to  (9). 

12.  After  reaching  the  last  downstream  station,  solve  the  pressure  equation 
(61)  by  tridiagnonal  matrix  algorithm.  Several  iterations  from  down¬ 
stream  to  upstream  are  employed  to  update  the  three-dimensional  ellip¬ 
tic  pressure  field. 

13.  Repeat  steps  (3)  to  (11)  for  several  sweeps  until  both  the  pressure  and 
velocity  fields  have  converged  within  a  specified  tolerance. 

14.  Return  to  step  (2)  for  the  next  time  step. 
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15.  Stop  if  the  steady-state  solution  is  achieved,  or  if  time  exceeds  the 
maximum  time  period  assigned.  For  steady-flow  calculations,  one  may 
relax  the  convergence  criterion  in  step  (1?)  and  use  a  larger  time 
increment  for  the  intermediate  solutions. 

In  the  present  study,  only  one  sweep  was  used  in  step  (12)  for  each  time 
step.  Also,  as  noted  earlier,  instead  of  specifying  the  initial  conditions 
for  velocity  and  turbulence  profiles  everywhere  at  t  =  0,  only  the  profiles  at 
the  first  station  5  =  1  are  specified,  and  the  downstream  profiles  are  taken 
from  the  immediately  upstream  station  (i.e.,  4> p  =  ^  at  t  =  0)  during  the 
first  global  sweep. 

IV.  SOME  NUMERICAL  ASPECTS  OF  THE  METHOD 

It  was  noted  in  Section  II  that  many  components  of  the  present  calcula¬ 
tion  method  have  been  tested  by  applying  it  to  study  the  flow  past  geometri¬ 
cally  simpler  bodies.  In  this  report,  therefore,  we  will  be  specifically 
concerned  with  the  flow  around  ship  hulls. 

Before  discussing  the  selection  of  specific  hull  forms  as  test  cases  and 
presenting  the  corresponding  results,  it  is  useful  to  evaluate  and  document 
certain  numerical  aspects  of  the  method.  These  are  discussed  in  relation  to 
the  calculations  performed  for  one  of  the  test  cases,  namely  the  SSPA  Liner, 
the  physical  aspects  of  which  will  be  considered  in  the  next  section.  For  the 
purposes  of  the  present  section  it  is  simply  a  representative  hull  form  which 
is  used  to  examine  the  convergence  properties  of  the  algorithms,  dependence  of 
the  solutions  on  the  grid  and  size  of  the  numerical  solution  domain,  and  the 
computer  times  required  to  obtain  acceptably  accurate  solutions. 

IV. 1  CONVERGENCE 

Two  important  measures  of  the  performance  of  an  iterative  numerical 
method  are  the  number  of  iterations  required  to  obtain  a  converged  solution 
and  the  influence  of  the  grid.  Here,  we  shall  examine  first  the  convergence 
history  of  several  representative  quantities  for  a  typical  calculation  with  a 
(50  x  30  x  15)  grid  in  the  (x,r,0)  directions.  In  this  particular  applica¬ 
tion,  as  well  as  in  all  others,  the  calculations  were  performed  for  160  time 
steps  (or  sweeps)  to  assure  full  convergence  of  the  solutions  although  much 
fewer  sweeps  are  actually  required  to  obtain  an  acceptable  solution. 
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Figure  7  shows  the  convergence  of  the  distributions  of  pressure,  Cp  (= 
2p),  and  friction  velocity!^,  along  the  keel  and  waterline  coordinates  and 
their  extensions  into  the  wake,  and  wake  centerline  velocity  Uc.  It  is  seen 
that  all  these  quantities  converge  monotonically  in  less  than  60  time  steps  or 
global  sweeps.  Recall  that  the  solutions  were  started  with  a  constant  ambient 
pressure  (p=0)  throughout  the  solution  domain.  The  calculated  pressure  and 
velocity  fields  clearly  capture  all  the  important  features  of  the  final  con¬ 
verged  solutions  in  less  than  20  sweeps.  This  is  particularly  encouraging 
because  it  eVninates  the  need  for  generating  an  initial  pressure  field  from  a 
potential-flow  calculation. 

Othe'  tests  of  convergence  of  the  solutions  made  on  the  mass  source  terms 
appearing  in  the  pressure  and  pressure-correction  equations  also  showed  mono¬ 
tonic  convergence  similar  to  that  depicted  for  the  flow  parameters  shown  in 
the  above  figures. 

IV. 2  GRID  DEPENDENCE 

Calculations  were  performed  with  four  different  (x,r,e)  grids  in  the 
solution  domain:  (0.3  <  x  <  4.5,  rs  <  r  <  1.0,  0  <  0  <  u/2) ,  to  examine  the 
sensitivity  of  the  numerical  solutions  to  grid  refinement.  Some  information 
on  these  calculations  is  given  in  Table  1,  and  typical  views  of  the  four  grids 
are  shown  in  Figure  8.  In  the  stern  and  near  wake  regions,  the  finest  grid, 
(74  x  30  x  21),  has  nearly  twice  as  many  nodes  in  each  direction  compared  with 
^e  coarsest  one,  (41  x  15  x  12).  The  coarsest  grid  calculation  converged  in 
only  15  iterations  and  required  only  about  40  cpu  seconds  on  the  CRAY 
XMP/48.  On  the  other  hand,  the  finest  grid  required  about  30  cpu  minutes  to 
achieve  the  same  level  of  convergence.  This  is  due  to  the  significant 
increase  of  the  number  of  grid  nodes  as  well  as  the  number  of  global  sweeps 
(110)  needed  for  convergence. 
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Table  1.  Sunmary  of  Grid  Dependence  Tests 


CASE  ==> 

I 

II 

III 

IV 

Grid  points 

in  x-direction 

41 

74 

50 

74 

in  r-direction 

15 

15 

30 

30 

in  a-direction 

12 

12 

15 

21 

Total  nodes 

7,380 

13,320 

22,500 

46,620 

Time  Step,  t 

1.0 

0.5 

0.3 

0.5 

Memory,  10^ 

0.32 

0.51 

0.84 

1.63 

cpu,  sec/iteration 

2.6 

5.3 

6.0 

16.6 

Iterations  for 

15 

60 

55 

no 

convergence 

Total  cpu  (secs) 

40 

320 

330 

1830 

Figure  9  shows  the  results  obtained  with  the  four  grid  arrangements  with 
regard  to  the  flow  parameters  considered  before.  It  is  quite  evident  that  the 
three  finer  grids  yield  nearly  identical  results.  On  the  other  hand,  the 
solution  with  the  coarsest  grid  shows  some  departure  from  the  other  three. 
This  is  presumably  due  to  an  inadequate  numerical  resolution  and  less  accurate 
specifications  of  boundary  conditions  at  the  hull  surface  and  along  the  wake 
centerline.  Nevertheless,  it  is  quite  encouraging  to  note  that  even  the 
coarsest  grid  calculation  is  able  to  capture  many  of  the  important  features  of 
the  flow  which  are  evident  from  the  data  discussed  later.  Consequently,  the 
coarse  grid  solutions  may  be  used  to  guide  the  selection  of  principal  param- 
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eters  in  a  design  process.  Refined  calculations  can  then  be  made  in  regions 
where  significant  variations  in  body  geometry  or  flow  features  occur. 

IV. 3  SOLUTION  DOMAIN 

In  some  earlier  studies  of  well  known  two-dimensional,  axi symmetric,  and 
simple  three-dimensional  bodies,  it  was  found  that,  with  the  simple  uniform 
stream  boundary  conditions  employed  here  the  numerical  solutions  were  quite 
sensitive  to  the  size  of  the  solution  domain.  Thp  most  critical  quantity  to 
be  examined  in  this  regard  is  the  pressure  distr  i..ution  since  its  influence 
penetrates  much  farther  into  the  inviscid-f low  region  compared  to  the  velocity 
field.  In  order  to  properly  capture  the  entire  zone  of  viscous-inviscid 
interaction,  the  outer  boundary  should  be  placed  at  a  distance  sufficiently 
far  away  from  the  body  such  that  the  uniform  flow  and  zero  pressure  conditions 
are  indeed  appropriate.  Similarly,  the  downstream  boundary  must  be  located  in 
the  far  wake  where  upstream  propagation  of  pressure  becomes  negligible. 

For  the  present  applications  to  ship  stern  and  wake  flows,  we  have  varied 
the  location  of  the  outer  and  downstream  boundaries  over  a  wide  range  to 
examine  the  influence  of  domain  size  on  the  solutions.  Calculations  were 
performed  for  four  different  combinations  of  downstream  and  outer  boundaries, 
namely 


0.3  <  x  <  23.1, 
0.3  <  x  <  4.53, 
0.3  <  x  <  1.95, 
0.3  <  x  <  1.34, 


rs  <  r  <  2.00 
rs  <  r  <  0.95 
rs  <  r  <  0.47 
rs  <  r  <  0.22 


with  a 

(57  x 

with  a 

(50  x 

with  a 

(45  x 

with  a 

(41  x 

35  x  15)  grid, 
31  x  15)  grid, 
27  x  15)  grid, 
22  x  15)  grid. 


The  grids  for  the  three  smaller  domains  were  obtained  by  simply  deleting  an 
appropriate  number  of  outer  and  downstream  grid  lines  from  the  grid  generated 
for  the  largest  solution  domain.  Therefore,  the  four  solutions  correspond  to 
essentially  the  same  grid  distributions.  It  is  seen  from  Figure  10  that  the 
solutions  with  the  two  larger  domains  are  essentially  the  same.  It  is  also 
clear  that  the  smallest  domain  is  too  small  to  correctly  predict  the  pressure 
distribution  over  the  hull.  These  calculations  suggest  that  domain  dependency 
can  be  eliminated  by  choosing  an  outer  boundary  which  is  farther  than  about  a 
half  ship  length  from  the  axis,  and  the  downstream  boundary  at  a  similar 


51 


distance  from  the  stern.  For  all  remaining  calculations  we  have  chosen  the 
domain  0.3  <  x  <  4.5,  rs  <  r  <  1.0  ,  which  is  more  than  adequate  to  capture 
the  viscous-inviscid  interaction. 


IV. 4  COMMUTING  TIMES 

The  initial  development  of  the  present  computer  code  and  some  preliminary 
calculations  were  performed  on  a  Prime  9950  minicomputer  of  the  University  of 
Iowa.  Quite  encouraging  results  were  obtained  even  though  the  finest  grid 
that  could  be  tested  was  just  (30  x  19  x  14).  However,  H  also  became  clear 
that  a  computer  with  greater  speed  and  larger  memory  was  needed  in  order  to 
realize  the  full  potential  of  the  method.  Therefore,  subsequent  studies  in 
grid  refinement,  domain  dependence,  and  convergence  were  performed  on  the  CRAY 
XMP/24  supercomputer  of  the  Naval  Research  Laboratory  (NRL)  and  the  CRAY 
XMP/48  machine  of  the  National  Center  for  Supercomputing  Applications  (NCSA) 
at  the  University  of  Illinois. 

In  order  to  fully  utilize  the  vectorization  capabilities  of  the  CRAY 
supercomputers,  several  major  revisions  had  to  be  made  in  the  original  code. 
These  included  vectorization  of  the  tridiagonal  algorithm  and  the  subroutines 
for  the  calculation  of  the  finite-analytic  coefficients.  These  code  optimiza¬ 
tions  resulted  in  a  70-percent  saving  in  cpu  time  compared  to  the  correspond¬ 
ing  scalar  calculations.  The  optimized  code  on  the  CRAY  supercomputer  runs 
about  150  times  faster  than  the  origina1  one  on  a  Prime  9950.  A  typical  ship 
stern  and  wake  flow  calculation  with  a  (50  x  30  x  15)  grid  now  takes  about  5- 
10  minutes  of  cpu  time  to  obtain  fully  converged  solutions.  It  is  estimated 
that  the  same  calculations  would  require  almost  15  hours  of  cpu  time  on  the 
Prime  9950. 


V.  TEST  CASES  AND  COMPARISONS  WITH  EXPERIMENTS 

One  of  the  major  problems  of  assessing  the  performance  of  a  complex 
computational  method  which  is  composed  of  many  numerical  and  physical  models 
and  approximations  is  to  find  sufficiently  detailed  and  varied  test  cases  for 
which  reliable  information  is  available  from  other,  independent  sources.  For 
turbulent  flows,  this  has  invariably  involved  recourse  to  experimental  infor¬ 
mation.  Ideally,  it  is  desirable  to  have  available  experimental  data  sets 
which  are  (a)  reliable  with  respect  to  accuracy,  (b)  complete  enough  to  pro- 
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vide  the  necessary  initial  and  boundary  conditions  required  to  carry  out  a 
meaningful  calculation,  and  (c)  contain  sufficient  additional  information 
which  could  be  used  to  evaluate  the  performance  of  the  calculation  method  and 
its  components.  These  requirements  were  brought  sharply  into  focus  at  the 
1980-81  Stanford  Conferences  on  Complex  Turbulent  Shear  Flows  at  which 
attempts  were  made  to  identify  such  data  sets  for  the  evaluation  of  turbulence 
models.  In  fact,  it  became  quite  evident  at  that  time  that  few  such  data  sets 
existed.  This  was  particularly  the  case  for  three-dimensional  flows. 

The  situation  with  regard  to  the  present  topic  of  ship  stern  and  wake 
flow  is  not  very  different.  In  spite  of  the  numerous  experimental  studies 
which  were  mentioned  in  Section  1 1.2,  there  is  not  a  single  data  set  which 
satisfies  all  of  the  requirements  noted  above.  However,  this  does  not  mean 
that  the  available  information  cannot  be  used  to  provide  guidance  in  the 
development  and  evaluation  of  calculation  methods.  If  the  limitations  are 
recognized,  the  data  can  still  be  used  to  evaluate  certain  aspects  of  the 
overall  computation  procedure.  This  also  implies  that  it  is  desirable  to  make 
comparisons  with  many  sets  of  data,  obtained  on  different  hulls,  by  different 
techniques,  and  in  different  facilities.  In  view  of  this,  we  have  selected 
five  cases  for  which  reasonably  extensive  data  are  available.  In  what  fol¬ 
lows,  we  shall  briefly  present  the  rational  for  selecting  each  of  them. 

The  Wigley  parabolic  hull  was  selected  as  the  first  test  case  because  it 
has  been  used  in  many  previous  studies  of  wave  and  viscous  resistance  in  ship 
hydrodynamics  and  because  it  has  a  rather  simple  stern  shape.  The  latter 
feature  is  attractive  in  the  evaluation  of  calculation  methods  with  respect  to 
the  prediction  of  the  very  near  wake  and  the  subsequent  evolution  of  the 
three-dimensional  wake. 

As  noted  earlier,  the  Wigley  hull  was  one  of  four  scapes  selected  for  the 
ITTC  Cooperative  Experimental  Program  which  sought  to  establish  a  comprehen¬ 
sive  data  base,  and  one  of  two  on  which  extensive  measurements  were 
reported.  However,  as  noted  in  ITTC  (1987),  many  of  these  experiments  were 
restricted  only  to  force  components.  Among  the  experiments  concerned  with  the 
details  of  the  viscous  flow  over  the  stern  and  in  the  wake,  the  most  complete 
are  those  of  Sarda  (1986)  on  a  double  model  in  a  wind  tunnel.  In  fact,  these 
experiments  were  carried  out  specifically  to  study  the  stern  and  wake  flow, 
and  document  the  evolution  of  the  wake  from  the  thin  boundary  layer  on  the 


53 


hull.  The  measurements  included  surface  pressure  distribution,  mean  velocity 
field,  and  all  components  of  the  Reynolds-stress  tensor.  The  comparisons 
between  these  experiments  and  the  present  computations  are  greatly  facilitated 
by  the  fact  that  the  data  were  collected  in  a  Cartesian  coordinate  system  and 
are  available  on  tape. 

The  second  test  case  is  the  SSPA  720  Cargo  Liner.  The  boundary  layer  on 
a  double  model  of  this  hull  form  was  measured  in  a  wind  tunnel  by  Larsson 
(1974).  These  are  among  the  most  well  documented  mean-flow  measurements  in 
three-dimensional  boundary  layers,  and  were  carefully  examined  and  recompiled 
in  a  convenient  form  in  preparation  for  the  SSPA-ITTC  Workshop  on  Ship  Bound¬ 
ary  Layers  (Larsson,  1981).  These  data  were  supplemented  by  corresponding 
turbulence  measurements  by  Lofdahl  (1982)  and  Lofdahl  and  Larsson  (1984). 
Because  these  experiments  were  conducted  to  study  the  hull  boundary  layer, 
rather  than  the  stern  and  wake  flow,  the  measurements  were  made  in  a  coordin¬ 
ate  system  suitable  for  three-dimensional  boundary  layers.  As  discussed  by 
Chen  and  Patel  (1985b),  this  makes  it  difficult  to  compare  the  data  with 
calculations  which  use  generalized  coordinates  because  the  errors  involved  in 
the  interpolations  are  not  insignificant.  Nevertheless,  the  data  of  Larsson 
and  Lofdahl  are  most  useful  in  scrutinizing  some  critical  aspects  of  the 
calculation  method,  and  also  in  gaging  the  advances  that  have  been  made  over 
the  methods  which  led  to  the  results  shown  in  Figure  1(a). 

The  HSVA  Tanker  is  the  third  test  case.  Mean  flow  measurements  in  the 
boundary  layer  on  a  double  model  of  this  hull  were  first  made  in  a  wind  tunnel 
by  Hoffman  (1976)  and  other  measurements  on  the  same  model  have  been  reported 
since  then.  For  example,  quite  detailed  mean-velocity  field  data  were 
obtained  by  Wieghardt  and  Kux  (1980)  and  Wieghardt  (1982,  1983),  and  some 
turbulence  measurements  were  reported  by  Knaack  (1984),  and  Knaack,  Kux  and 
Wieghardt  (1985).  Recent  communications  with  these  authors  indicate  that 
other  experiments  on  this  hull  are  still  continuing. 

For  the  comparisons  made  in  this  report,  we  have  used  the  tabulated  stern 
and  wake  flow  data  supplied  on  tape  by  Dr.  Kux.  Comparison  of  the  present 
calculations  with  these  data  was  greatly  facilitated  by  the  fact  that  the 
calculated  results  could  be  readily  interpolated  into  the  Cartesian  coordin¬ 
ates  in  which  the  measurements  were  made.  The  earlier  boundary-layer  measure¬ 
ments  of  Hoffman  were  used  only  as  a  guide  to  determine  appropriate  initial 
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conditions.  Recall  that  these  measurements  had  led  to  the  selection  of  this 
hull  as  the  second  test  case  at  the  SSPA-ITTC  Workshop.  Comparison  of  the 
present  results  with  those  presented  at  the  Workshop  (Figure  lb)  would  again 
Indicate  the  progress  that  has  been  made  in  the  prediction  of  stern  and  wake 
flows.  Futhermore,  the  HSVA  Tanker  is  of  fuller  form  than  the  SSPA  hull  and 
therefore  provides  a  more  severe  test  of  the  calculation  method. 

The  next  test  case  is  somewhat  similar  to  the  previous  one  insofar  as  it 
also  involves  a  full  form  hull  and  mean-flow  measurements  on  a  double  model  in 
a  wind  tunnel.  This  is  the  SRI 07  Ore  Carrier  which  has  been  the  subject  of 
many  collaborative  experiments  in  Japan  (see,  for  example,  Nagamatsu,  1981; 
Okajima,  Toda,  and  Suzuki,  1985).  Although  these  data  are  restricted  to  the 
stern  and  wake  region,  and  therefore  there  is  some  uncertainty  concerning  the 
proper  initial  conditions,  they  were  selected  because  they  enable  the  eval¬ 
uation  of  the  calculation  method  against  somewhat  similar  data  sets  obtained 
in  quite  different  facilities.  Here  again,  the  use  of  Cartesian  coordinates 
in  the  experiments  made  it  convenient  to  carry  out  detailed  comparisons. 

The  fifth  and  final  case  considered  is  the  well  known  Series  60,  Cb  = 
0.60  hull.  This  is  the  second  hull  for  which  the  ITTC  Cooperative  Experimen¬ 
tal  Program  led  to  the  collection  of  extensive  towing- tank  data.  As  in  the 
case  of  the  Wig  ley  hull,  however,  most  of  these  experiments  were  concerned 
with  force  components  and  none  documented  the  stern  and  wake  flow  in  as  much 
detail  as  that  on  the  other  hulls  mentioned  above.  Fortunately,  quite  exten¬ 
sive  data  for  this  hull  have  recently  become  available  from  an  independent 
study. 

Here,  we  shall  use  the  data  obtained  in  the  course  of  a  joint  research 
program  on  propeller-hull  interaction  between  Osaka  University  and  the 
University  of  Iowa.  These  experiments  included  measurements  with  and  without 
an  operating  propeller  and  will  be  described  in  a  forthcoming  report  by  Toda, 
et  al.,  (1988).  For  the  present  purposes,  we  will  consider  only  the  measure¬ 
ments  without  the  propeller.  These  experiments  on  the  Series  60  are  particu¬ 
larly  suitable  for  the  present  study  for  two  reasons.  First,  they  were 
designed  specifically  with  the  goal  of  obtaining  data  for  the  validation  of 
modern  computational  methods,  and  therefore  the  measurement  coordinates  and 
locations  were  selected  to  facilitate  direct  comparisons  with  predictions. 
Secondly,  the  measurements  were  made  in  a  towing  tank  on  a  large  model  at  a 
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relatively  low  Froude  number,  and  therefore  comparisons  with  double-body 
calculations  give  some  information  on  the  influence  of  the  free  surface  in  the 
absence  of  a  strong  wave  system. 

Taken  together,  the  five  hull  forms  selected  here  provide  a  rather  com¬ 
prehensive  set  of  test  cases  for  the  validation  of  the  computational  method. 
Included  are  simple  and  complex  geometries  with  different  bow  and  stern 
shapes,  and  different  fullness,  data  obtained  on  double  models  in  different 
wind  tunnels  and  those  obtained  with  a  conventional  model  in  a  towing  tank, 
and  finally,  measurements  with  instrumentation  of  varying  degrees  of  sophisti¬ 
cation.  In  the  following  sections,  we  shall  point  out  further  features  of 
each  of  the  data  sets  which  have  a  bearing  on  the  evaluation  of  the  perfor¬ 
mance  of  calculation  procedures,  in  general,  and  the  present  method,  in  parti¬ 
cular.  In  the  calculations  presented  below,  all  attempts  have  been  made  to 
carefully  match  the  conditions  prevailing  in  the  experiment  but  lack  of  some 
critical  information  on  the  experiment,  on  the  one  hand,  and  limitations  of 
the  calculation  method,  on  the  other,  need  to  be  borne  in  mind  in  drawing 
general  conclusions  from  the  comparisons. 

Before  presenting  the  resits  For  the  individual  cases,  it  is  useful  to 
first  describe  some  features  that  are  common  to  all  of  them.  Among  these  are 
the  notation  used,  and  the  procedure  for  the  generation  of  initial  and  bound¬ 
ary  conditions. 

All  coordinates  and  geometrical  parameters  have  been  rendered  dimension¬ 
less  using  the  ship  length  L  as  the  characteristic  scale,  the  coordinate  x. 
being  measured  from  the  bow,  rather  than  from  midships.  The  velocity  compon¬ 
ents  are  nondimensional ized  by  a  constant  reference  velocity,  U0,  which  is 
either  the  velocity  of  the  model  or  the  velocity  in  the  tunnel  working  sec¬ 
tion.  Finally,  pressures  and  stresses  are  made  dimensionless  using  the  same 
reference  velocity,  and  ambient  density  and  viscosity.  All  other  quantities 
will  be  defined  when  they  are  first  introduced.  It  should  also  be  noted  that 
the  dimensionless  coordinates  (x,y,z)  are  used  interchangeably  with  upper  case 
notations  (X,Y,Z)  in  the  figures. 

The  initial  and  boundary  conditions  required  by  the  present  method  were 
discussed  in  general  terms  in  Section  1 1 1. 2  (d).  To  obtain  solutions  for  a 
particular  shape,  it  is  necessary  to  prescibe  the  conditions  at  the  upstream 
and  outer  surfaces  of  the  solution  domain.  Here,  we  have  employed  uniform- 
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stream  conditions  at  the  outer  boundary  (i.e.,  U  =  1,  W  =  0,  and  p  =  0)  and, 
therefore,  the  solutions  correspond  to  the  flow  around  double  bodies  in  an 
unrestricted  stream.  In  other  words,  wind-tunnel  or  towing-tank  blockage 
effects  are  not  taken  into  account. 

At  the  upstream  section,  x  =  0.3,  it  is  necessary  to  prescribe  the  dis¬ 
tributions  of  (U,V,W,k,e)  at  all  grid  points,  within  the  thin  boundary  layer 
as  well  as  in  the  inviscid  flow  outside.  As  explained  earlier,  it  is  not 
necessary  to  prescribe  the  pressure  because  it  is  determined  implicitly  from 
the  equation  of  continuity.  The  specification  of  the  five  quantities,  for 
which  there  are  corresponding  transport  equations,  requires  some  care  because 
they  must  properly  reflect  the  upstream  history  of  the  flow.  In  all  cases  we 
prescribe  a  girthwise  distribution  of  the  boundary-layer  thickness  6,  the 
friction  coefficient  Cf  =  2  u  ,  and  the  velocity  at  the  edge  of  the  boundary 
layer  U^.  These  are  estimated  from  either  previous  boundary- layer  and  invis- 
cid-flow  calculations,  or  guessed  and  then  adjusted  to  obtain  agreement  with 
data  at  the  first  measurement  station  downstream  from  the  initial  section. 
These  quantities  are  used,  together  with  the  law  of  the  wall  and  the  law  of 
the  wake,  to  generate  the  profiles  of  the  longitudinal  velocity  U  inside  the 
boundary  layer,  and  the  reduction  from  U6  to  unity  in  the  inviscid  flow  is 
assumed  to  take  place  as  r“2.  In  the  first  instance,  V  and  W  components  are 
set  to  zero,  and  k  and  e  are  obtained  from  correlations  for  a  flat-plate 
boundary  layer.  As  the  solution  progresses,  however,  the  values  of  V,  W,  k 
and  e  are  updated  by  scaling  those  calculated  at  the  first  downstream  sec¬ 
tion.  This  process  is  continued  only  for  the  first  20  global  sweeps,  and  then 
the  initial  profiles  are  fixed.  This  rather  intricate  procedure  for  the 
generation  of  the  initial  conditions  does  not  affect  the  principal  quantity, 
i.e.,  the  axial  velocity  profile.  However,  it  ensures  that  the  subsequent 
solution  is  carried  out  with  initial  profiles  of  transverse  velocity  compon¬ 
ents  and  turbulence  parameters  which  are  compatible  with  the  governing  equa¬ 
tions.  The  procedure  is  quite  automatic  and  results  in  solutions  which  are 
not  sensitive  to  initial  conditions. 


57 


V.l  WIGLEY  PARABOLIC  HULL 


(a)  Available  Experimental  Information 

The  experiments  of  Sarda  (1986)  were  carried  out  on  a  10  ft  (3.048  m) 
long  double  model  suspended  by  cables  in  a  5  ft  (1.52  m)  octagonal-section, 
return-circuit,  closed  wind  tunnel  in  which  the  freestream  turbulence  level 
was  less  than  0.5  percent.  The  measurements  were  made  at  a  Reynolds  number  Re 
=  UqL/v  =  4.5  x  10®.  Although  measurements  could  be  made  only  upto  a  distance 
of  8  ft  (2.44  m,  or  0.8  L)  downstream  of  the  stern,  these  experiments  document 
the  wake  considerably  farther  than  the  other  experiments  discussed 
subsequently.  The  boundary  layer  on  the  hull  was  tripped  by  means  of  studs  of 
standard  design  at  x  =  0.05. 

The  available  data  include  the  pressure  distribution  on  the  hull  measured 
by  surface  pressure  taps,  mean-velocity  components  measured  by  a  five-hole 
pitot  probe,  and  components  of  the  mean  velocity,  the  Reynolds-stress  tensor, 
and  triple  products  of  fluctuating  velocity  components  measured  by  means  of  a 
three-sensor  hotwire  probe.  The  velocity  and  turbulence  measurements  were 
made  at  several  longitudinal  sections  in  the  range  0.5  <  x  <  1.8,  i.e.,  from 
midships  to  0.8  ship  lengths  downstream  of  the  ste'-n. 

(b)  Body  Shape  and  Numerical  Grid 

The  coordinates  of  the  Wigley  hull  are  given  by 

2Y/B  =  [1  -  (2X/L  -  l)2]  [1  -  2Z/D)2]  (79) 

where  ()T,7,Z)  are  Cartesian  coordinates  with  the  origin  at  the  bow,  and  the 
ship  length  (L):  beam  (B):  draft  (D)  ratio  is  10:1:0.625. 

The  solution  domain  extends  from  a  section  0.3  ship  lengths  from  the  bow 
to  3.524  ship  lengths  downstream  of  the  stern,  and  from  the  hull  to  a  cylin¬ 
drical  surface  located  one  ship  length  from  the  ship  longitudinal  axis.  In 
terms  of  the  coordinate  variables,  the  solution  domain  is  defined  by  0.3  <  x 
<  4.524,  rs  <  r  <  1.0  .  This  domain  is  covered  by  50  x  30  x  15  grid  points  in 
the  axial,  radial,  and  circumferential  directions,  respectively,  giving  a 
total  of  22,500  grid  nodes. 
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Some  views  of  the  numerically-generated  grid  are  shown  in  Figure  11.  The 
grid-control  functions  f"1  used  to  obtain  this  grid  were  given  in  equation 
(41).  The  parameters  appearing  therein  are  as  follows: 

0.6  if  A  <  0 

0.2  if  A  >  0 

d  =  0.2 

The  following  are  some  of  the  noteworthy  features  of  this  grid.  First,  the 
coordinate  surface  n  =  2,  just  off  the  hull  in  the  radial  direction,  is 
arranged  in  such  a  way  that  at  all  locations  it  lies  in  the  range  90  <  y+  < 
250,  y+ 'being  the  normal  distance  from  the  hull.  Recall  that  such  a  con¬ 
straint  is  necessary  for  the  wall-functions  approach  to  be  applicable.  Other 
investigators  who  also  employ  wall  functions  have  placed  the  first  near-wall 
grid  surface  much  closer  to  the  hull  than  the  one  used  here.  While  such 
solutions  can  be  also  obtained  with  the  present  method,  we  believe  that  they 
cannot  be  trusted  because  the  underlying  assumptions  of  the  law  of  the  wall 
and  energy  equilibrium  are  no  longer  valid.  In  other  words,  grids  finer  than 
the  ones  used  here  can  be  justified  only  if  the  wall-function  approach  is 
abandoned  in  favor  of  direct  solutions  of  the  governing  equations  all  the  way 
upto  the  wall.  Secondly,  we  note  that  the  grid  is  concentrated  near  the  stern 
and  in  the  near-wake  to  resolve  the  large  gradients  which  occur  at  the  trail¬ 
ing  edge.  Third,  the  grid  spacing  is  increased  very  rapidly  in  the  radial  and 
longitudinal  directions.  In  fact,  the  final  step  size  in  the  radial  and  axial 
directions  is  of  the  order  of  0.2  and  0.8,  respectively.  The  use  of  such 
large  steps  while  maintaining  solution  accuracy  is  made  possible  by  the  fin¬ 
ite-analytic  discretization  of  the  transport  equations. 

(c)  Description  of  Results  and  Comparisons  with  Experiments 

Here,  as  in  all  subsequent  cases,  we  shall  first  examine  the  distribu¬ 
tions  of  pressure  and  wall  shear  stress  on  the  hull  surface,  and  then  the  de¬ 
tails  of  the  flow  field.  Figure  12  shows  the  distribution  of  pressure  coeffi¬ 
cient  CD  (=2p)  and  friction  velocity  U  along  the  keel  and  the  waterline 

planes  of  synmetry.  Girthwise  distributions  of  Cn  and  the  magnitude  of  the 

2  ^ 

friction  coefficient  Cf  (=  2IJT)  at  a  few  transverse  sections  are  shown  in 
Figure  13.  The  data  of  Sarda  (1986)  and  some  unpublished  results  from  Watmuff 
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and  Joubert  (1985)  are  also  shown  in  Figures  12  and  13.  In  the  case  of 
Watmuff  and  Joubert  the  pressure  datum  was  adjusted  to  match  the  potential 
flow  value  at  midships.  The  direction  of  the  wall  shear  stress  can  be  seen 
from  the  limiting  or  wall  streamlines  plotted  in  Figure  14.  This  last  figure 
was  obtained  from  the  calculated  wall  shear  stress  vectors  over  the  hull, 
starting  from  the  initial  section  at  x  =  0.3. 

In  Figure  12  the  experimental  data  are  restricted  to  the  waterline 
because  there  are  no  pressure  taps  along  the  sharp  keel.  The  potential-flow 
pressure  distributions  were  obtained  at  the  David  Taylor  Research  Center  using 
the  XYZ  computer  program.  It  is  seen  that  the  present  calculations  are  in 
good  agreement  with  the  data  of  Watmuff  and  Joubert,  but  both  are  considerably 
higher  than  the  measurements  of  Sarda  near  the  stern,  x  >  0.94.  The  discrep¬ 
ancy  between  the  two  sets  of  measurements  is  somewhat  surprising  because  both 
experiments  were  performed  under  similar  conditions.  Although  the  reason  for 
the  observed  differences  is  not  entirely  clear,  Sarda  noted  that  the  geometry 
of  his  model  was  not  precisely  the  same  as  that  of  the  mathematical  form,  his 
model  being  thicker  in  the  bow  and  stern  regions.  Potential-flow  calculations 
for  the  actual  experimental  shape  is  obviously  needed  to  determine  if  the 
differences  in  geometry  are  indeed  responsible  for  the  differences  in  the 
data.  Despite  these  uncertainties,  however,  we  note  that  the  present  results 
are  in  good  agreement  with  the  potential-flow  solutions  over  much  of  the  hull 
where  the  boundary  layer  is  thin.  The  level  of  disagreement  between  the  data 
of  Sarda  and  the  calculations  over  the  midbody  is  also  consistent  with  the 
expected  tunnel-blockage  effect  in  the  experiments.  In  the  stern  region,  the 
calculated  pressure  gradients  are  smaller  than  those  in  potential  flow  as  is 
to  be  expected  from  the  displacement  effect  of  the  viscous  flow. 

The  calculated  friction  velocities  along  the  waterline  follow  the  trends 
shown  by  the  data,  which  were  obtained  from  Clauser  plots  of  measured  velocity 
profiles,  although  the  measurements  are  somewhat  lower  than  the  predicted 
values.  It  is  seen  that  the  friction  velocity  decreases  along  the  waterline 
due  to  the  continuous  thickening  of  the  boundary  layer.  On  the  other  hand, 
the  calculated  friction  velocity  increases  along  the  keel  in  spite  of  the 
adverse  pressure  gradient  that  is  present.  This  is  due  to  a  reduction  in  the 
boundary  layer  thickness  produced  by  flow  divergence  out  of  the  keel  plane. 
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From  Figure  13  it  is  seen  that  the  calculated  girthwise  distributions  of 
pressure  are  in  agreement  with  the  two  sets  of  measurements  upto  x  =  0.95. 
The  marked  disagreement  between  the  calculations  and  the  data  of  Watmuff  ar.d 
Joubert  at  the  next  section,  x  =  0.962,  is  most  likely  due  to  interference 
from  the  model  support  in  the  experiment.  Further  downstream,  the  calcula¬ 
tions  agree  with  the  data  of  Watmuff  and  Joubert  but,  as  noted  already,  lie 
above  the  measurements  of  Sarda. 

The  calculations  and  experiments  indicate  that  in  the  region  0.5  <  x  < 
0.9  the  pressure  decreases  in  the  girthwise  direction  from  the  keel  towards 
the  waterline.  This  is  in  agreement  with  potential-flow  theory,  as  would  be 
expected  in  regions  where  the  boundary  layer  is  thin.  Although  the  girthwise 
gradient  is  small,  it  is  responsible  for  the  secondary  motion  from  the  keel 
towards  the  waterline.  Around  x  =  0.8,  the  pressure  shows  an  increase  towards 
the  waterline  in  a  region  whose  girthwise  extent  increases  downstream.  Even¬ 
tually,  for  x  >  0.9,  there  is  an  increase  in  pressure  all  the  way  from  the 
keel  to  the  waterline.  This  reversal  of  trend  from  the  predictions  of  invis- 
cid  theory  is  associated  with  the  changes  in  the  girthwise  distribution  of 
boundary  layer  thickness  over  the  stern. 

The  calculated  and  measured  wall  shear-stress  coefficients  shown  in 
Figure  13  indicate  similar  general  trends  but  there  exist  systematic  differ¬ 
ences,  particularly  at  the  most  upstream  station,  where  the  calculated  values 
are  higher  than  the  measurements.  Although  somewhat  better  agreement  could  be 
secured  in  this  respect  by  adjustments  in  the  initial  conditions  at  x  =  0.3, 
this  was  not  pursued  further  partly  because  the  agreement  between  the  calcula¬ 
tions  and  the  data  improves  downstream  (x  >  0.92).  The  decrease  ^n  wall  shear 
stress  magnitude  from  the  keel  towards  the  waterline  is  consistent  with  the 
thinning  of  the  boundary  layer  along  the  keel  and  a  thickening  along  the 
waterline  brought  about  by  the  crossflow.  The  general  direction  of  the  cross- 
flow  and  the  changes  in  the  direction  of  the  wall  shear  stress  are  best  seen 
from  the  wall  streamlines  of  Figure  14. 

An  overview  of  the  flow  pattern  is  provided  by  Figure  15  which  shows  the 
velocity  vectors  projected  onto  the  horizontal  waterplane  and  the  vertical 
centerplane.  These  plots  were  constructed  by  interpolation  of  the  numerical 
results.  Contours  of  the  axial  velocity  (U)  and  velocity  vectors  projected 
into  transverse  sections  (V  and  W  components)  are  shown  in  Figures  16(a)  and 
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(b)  for  some  representative  sections  at  midships,  over  the  stern,  in  the  near 
wake,  and  farther  downstream.  Here,  the  transverse  velocity  plots  also  indi¬ 
cate  the  grid  distribution.  We  recall  that  all  of  these  figures  show  informa¬ 
tion  close  to  the  hull  although  the  total  solution  domain  is  much  larger.  The 
measurements  of  Sarda  are  shown  in  Figures  17  (a)  and  (b)  for  comparison. 

Among  the  major  features  evident  from  Figures  15  and  16  are  the  follow¬ 
ing.  In  the  longitudinal  direction,  the  boundary  layer  grows  rapidly  along 
the  waterline  but  remains  thin  along  the  keel.  This  is  due  to  the  secondary 
motion  which  is  directed  from  the  keel  towards  the  waterline,  in  the  viscous 
flow  as  well  as  in  the  outer  inviscid  flow.  The  transport  of  the  boundary 
layer  fluid  towards  the  waterline  is  also  indicated  by  the  convergence  of  the 
wall  streamlines  in  Figure  14.  The  direction  of  the  transverse  motion  is  also 
consistent  with  that  expected  from  considerations  of  the  pressure  gradients 
associated  with  an  inviscid  flow.  From  potential-flow  pressure  distribution 
on  the  front  half  of  the  hull  we  would  expect  the  secondary  flow  in  the  bound¬ 
ary  layer  upto  midships  and  perhaps  for  some  distance  downstream  to  be 
directed  towards  the  keel,  resulting  in  a  thicker  boundary  layer  at  the 
keel.  These  features  were  observed  in  the  measurements  of  Sarda  at  x  =  0.5 
(see  Figure  17)  and  also  in  his  boundary- layer  calculations.  The  present 
calculations  at  x  =  0.5  indicate  that  although  a  thicker  boundary  layer  is 
observed  at  the  keel  at  this  section,  the  secondary  motion  is  directed  towards 
the  keel  only  in  a  small  layer  close  to  the  keel.  This  difference,  although 
small,  is  most  likely  due  to  the  use  of  simple  flat-plate  profiles  for  initial 
conditions  (at  x  =  0.3)  in  the  present  calculation.  While  the  use  of  a  bound¬ 
ary-layer  calculation  to  obtain  the  initial  conditions  would  reduce  this 
difference,  additional  difficulties  arise  in  handling  the  sharp  keel  within 
the  framework  of  boundary-layer  theory. 

Examination  of  the  velocity  field  over  the  stern  and  in  the  near  wake 
shows  no  dramatic  changes  in  the  flow  structure.  In  particular,  there  is  no 
evidence  of  a  stern  vortex,  although  there  is  obviously  longitudinal  vorticity 
associated  with  the  transverse  velocity  components.  The  transverse  components 
decay  rather  rapid^v  in  the  wake.  The  calculations  predict  a  quite  ^apid 
recovery  of  the  wake  to  what  appears  to  be  en  almost  axisymmetric  state  by  x  = 
1.5,  a  half  ship  length  from  the  stern.  This  is  at  variance  with  the  measure¬ 
ments  shown  in  Figure  17  which  indicate  a  rather  slow  development  of  the  wake. 
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The  calculated  distributions  of  the  velocity  components  (U,V,W),  pressure 
(p),  and  turbulent  kinetic  energy  (k)  are  compared  with  the  measurements  of 
Sarda  (1986)  in  Figure  18  at  several  sections  in  the  range  0.5  <  x  <  1.5.  The 
calculated  results  were  interpolated  to  determine  the  variation  of  the  various 
quantities  in  the  horizontal  (y)  direction  at  different  depths  (z  =  constant), 
z  =  0  being  the  waterplane.  Note  that  the  horizontal  axis  in  these  figures  is 
y*,  the  distance  from  the  hull  surface.  For  profiles  below  the  keel  (z  < 
-0.0625)  and  in  the  wake  (x  >  1),  of  course,  y*  =  y.  The  data  obtained  with  a 
five-hole  pitot  and  a  triple-sensor  hotwire  are  included  to  show  the  exper¬ 
imental  uncertainty,  particularly  with  respect  to  the  components  V  and  W.  In 
the  case  of  pressure,  the  original  data  which  were  obtained  with  the  pitot 
probe  were  referenced  to  a  pressure  in  the  wind  tunnel  some  distance  upstream 
of  the  model.  The  pressure  measured  in  the  wake  at  x  =  1.5  was  nearly  uniform 
but  higher  than  the  reference  pressure.  The  calculations,  on  the  other  hand, 
predicted  almost  zero  pressure  relative  to  the  pressure  in  the  uniform 
stream.  In  view  of  this,  a  constant  (=  -0.06)  reference-pressure  correction 
has  been  added  to  the  data  to  obtain  the  so-called  corrected  pressures  shown 
in  Figure  18.  In  other  words,  the  calculated  and  measured  pressures  have  been 
matched  in  the  wake  at  x  =  1.5.  Such  a  constant  correction  obviously  does  not 
alter  the  pressure  gradients  and  therefore  enables  us  to  compare  the  calcu¬ 
lated  and  measured  gradients.  On  the  other  hand,  it  is  not  a  correction  for 
tunnel  blockage  although  blockage  may  be  responsible  for  a  part  of  the 
observed  pressure  difference.  A  more  likely  possibility  is  an  error  in  probe 
calibration.  Figure  18  contains  a  great  deal  of  information  regarding  the 
details  of  the  flow  and  on  the  performance  of  the  calculation  method.  To  aid 
the  understanding  of  the  many  issues  involv.d,  we  shall  discuss  the  results 
starting  with  those  at  midships  and  progressing  downstream  into  the  wake. 

The  results  a'  midships  (Figure  18a)  indicate  that  the  boundary  layer  is 
thin  and  the  girthwise  variation  of  its  thickness  is  captured  adequately  by 
the  simple  initial  conditions  prescribed  at  x  =  0.3.  There  are  differences  in 
the  shapes  of  the  axial  velocity  profiles  at  some  depths,  and  in  the  shapes 
and  magnitudes  of  the  transverse  components  particularly  close  to  the  wall. 
In  spite  of  the  differences  between  the  pitot  and  hotwire  data,  it  is  clear 
that  the  generally  negative  values  of  V  and  W  measured  in  the  experiments  are 
not  captured  by  the  calculations.  We  have  already  commented  on  this  differ¬ 
ence.  It  appears  that  further  refinements  in  the  initial  conditions  could 


have  been  made  to  secure  better  agreement  at  this  first  measurement  section. 
In  spite  of  the  differences  in  the  mean  velocity  field,  we  see  that  the  cal¬ 
culated  pressure  is  in  reasonable  agreement  with  the  corrected  data.  The 
small  departures  between  the  two  close  to  the  wall  are  believed  to  be  due  to 
probe  interference.  Finally,  the  profiles  of  turbulent  kinetic  energy  are 
also  matched  rather  well  except  at  the  keel  (z  =  -0.0625). 

The  flow  over  the  stern  is  described  by  the  results  at  four  sections,  x  = 
0.90,  0.95.  0.967  and  0.983.  Here,  the  profiles  of  the  axial  velocity  indi¬ 
cate  that  the  girthwise  distribution  of  the  thickness  of  the  viscous  region  is 
predicted  rather  well,  including  the  near-disappearance  of  the  boundary  layer 
at  the  keel.  There  is,  however,  a  tendency  for  the  calculated  velocities  to 
be  larger  close  to  the  wall.  This  is  most  likely  associated  with  the  use  of 
the  wall  functions.  The  correction  for  pressure  gradients  in  the  law  of  the 
wall  are  presumably  not  sufficient  to  <~^mpletely  describe  the  wall  layer. 
With  respect  to  the  transverse  components  of  velocity,  the  calculations  pre¬ 
dict  the  general  trends  and  magnitudes  rather  well.  Unfortunately  the  scatter 
in  the  data  and  the  differences  between  the  pitot  and  hotwire  results  preclude 
more  definitive  conclusions.  The  predicted  distributions  of  pressure  are  in 
good  agreement  with  the  corrected  data  with  the  possible  exception  of  a  short 
fetch  close  to  the  stern.  The  profiles  of  turbulent  kinetic  energy  are  also 
predicted  quite  well  except  very  close  to  the  wall  where  the  measured  values 
are  generally  larger.  This  feature  is  presumaby  related  to  the  differences  in 
the  axial  velocity  profiles  noted  above  and  may  again  be  due  to  the  inadequacy 
of  the  wall  functions. 

The  results  at  the  next  three  stations,  namely,  x  =  1.002,  1.017  and 
1.050,  illustrate  the  rapid  changes  that  take  place  in  the  neighborhood  of  the 
centerplane  (y*  =  y  =  0)  of  the  extreme  near  wake.  First,  we  note  that  the 
calculated  results  dc  not  extend  to  the  centerplane  because  of  the  way  in 
which  the  plane-of-symmetry  boundary  conditions  are  enforced.  However,  this 
is  not  the  principal  reason  for  the  quite  major  differences  between  the  calcu¬ 
lations  and  measurements.  Previous  calculations  for  a  flat  plate  (Patel  and 
Chen,  1987)  and  for  bodies  of  revolution  (Chen  and  Patel, 1987b,  1988)  Indicate 
similar  tendencies,  i.e.,  an  overestimation  of  the  velocity  and  an  underesti¬ 
mation  of  the  turbulent  kinetic  energy  along  the  wake  centerline.  This  is 
attributed  to  the  fact  that  the  flow  in  the  sublayer  and  the  blending  zone  of 
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the  upstream  boundary  layer  has  not  been  considered  as  a  result  of  the  appli¬ 
cation  of  wall  functions  in  the  tubulence  model.  The  lower  measured  veloci¬ 
ties  and  higher  measured  turbulence  are  therefore  associated  with  the  physics 
of  the  flow  in  the  near-wall  region.  The  secondary  velocity  components  also 
indicate  systematic  differences,  presumably  due  to  the  same  reason.  In  spite 
of  the  discrepancy  in  the  inner  portion  of  the  very  near  wake,  it  is  encourag¬ 
ing  to  note  that  the  flow  in  the  outer  part  continues  to  be  predicted  satis¬ 
factorily  in  all  respects. 

The  subsequent  evolution  of  the  wake  is  depicted  by  the  results  at  the 
next  four  sections,  x  =  1.1,  1.2,  1.3  and  1.5.  It  is  clear  that  the  differ¬ 
ences  which  arose  in  the  inner  part  of  the  wake  at  the  stern  persist  for  quite 
large  distances,  particularly  with  regard  to  the  profiles  of  axial  velocity 
and  turbulent  kinetic  energy.  The  scatter  in  the  data  is  such  that  definitive 
conclusions  concerning  the  secondary  motion  cannot  be  reached  but  there  is 
some  evidence  to  suggest  that  the  calculations  predict  a  more  rapid  decay  of 
the  secondary  flow  than  is  observed  in  the  data.  This  was  also  seen  from 
comparisons  of  Figures  16  and  17.  Finally,  we  observe  that  the  pressure  at 
these  distances  is  nearly  uniform,  and  the  comparison  between  the  data  and  the 
calculations  for  x  =  1.5  show  that  a  constant  correction  on  the  measured 
pressures  is  justified. 

(d)  Summary 

Insofar  as  the  flow  over  the  hull  and  the  stern  is  concerned,  it  is  clear 
that  the  present  calculations  predict  all  essential  features  of  the  mean  flow 
with  reasonable  accuracy.  No  attempt  has  been  made  here  to  compare  the  pre¬ 
sent  results  with  those  obtained  earlier  by  solutions  of  the  boundary-layer 
equations  (e.g.  by  Sarda,  1986),  but  it  is  expected  that  the  differences  may 
be  significant  only  in  the  stern  region  because  there  is  only  a  weak  viscous- 
inviscid  interaction  on  this  rather  slender  hull.  With  regard  to  the  turbu¬ 
lence,  we  have  chosen  to  examine  only  the  turbulent  kinetic  energy,  a  quantity 
which  can  be  measured  with  a  greater  degree  of  precision  than  the  individual 
Reynolds-stress  components,  and  which  directly  enters  the  turbulence  model 
used  in  the  method.  The  results  could  of  course  be  analyzed  further  to  exam¬ 
ine  other  aspects  of  turbulence  modelling.  This  has  not  been  pursued  in  view 
of  the  following. 
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Perhaps  the  most  interesting  aspect  of  the  calculations  for  the  Wigley 
hull  is  the  very  large  and  marked  influence  of  using  the  wall  functions  in  the 
turbulence  model  on  the  prediction  of  the  wake.  The  experiments  indicate  that 
the  three  dimensionality  and  turbulence  generated  in  the  near-wall  region  of 
the  stern  flow  persist  for  quite  large  distances  in  the  wake.  This  is  presum¬ 
ably  because  the  pressure  gradients  in  the  stern  region  of  this  hull  are  too 
weak  to  significantly  affect  the  development  of  the  near  wake.  It  is  clear 
that  further  improvements  in  the  prediction  of  the  wake  would  require  replace¬ 
ment  of  the  wall  functions  with  a  turbulence  model  that  would  enable  accurate 
resolution  of  the  near-wall  layers. 

V.2  SSPA  CARGO  LINER 

(a)  Hull  Geometry  and  Available  Experimental  Information 

The  experiments  of  Larsson  (1974)  were  conducted  on  a  double  model  of  a 
cargo  liner  in  a  1.25  m  x  1.80  m  wind  tunnel,  with  a  freestream  turbulence 
level  of  about  0.2  percent.  The  model  was  supported  from  the  tunnel  floor  by 
two  struts.  The  tests  were  conducted  at  a  Reynolds  number  of  5  x  10  ,  based 
on  model  length  L  =  2  m.  The  hull  offsets  are  shown  in  Figure  19.  The  beam 
and  draft  of  the  model  are  0.283  m  and  0.118  m,  respectively,  and  the  block 
coefficient  is  0.675.  Larsson  made  extensive  potential-flow  calculations, 
with  and  without  tunnel  walls,  to  study  the  effect  of  wind-tunnel  blockage  on 
the  hull  pressure  distribution.  These  calculations  have  been  used  here  to 
correct  the  measured  pressures  for  blockage. 

As  this  was  an  investigation  of  the  hull  boundary  layer,  measurements  of 
velocity  profiles  were  made  in  streamline  coordinates,  i.e.,  the  two  velocity 
components  in  planes  parallel  to  the  hull  surface,  along  and  normal  to  the 
inviscid-f low  streamlines  on  the  hull,  were  measured  along  local  normals  to 
the  hull.  The  inviscid  streamlines  calculated  by  Larsson,  and  along  which  his 
measurements  were  made,  are  shown  in  Figure  19(c).  The  mean-velocity  measure¬ 
ments  were  made  by  hotwires  and  covered  the  range  0.20  <  x  <  0.95.  In  a  later 
investigation,  Lofdahl  (1982)  (see  also  Lofdahl  and  Larsson,  1984)  measured 
all  six  components  of  the  Reynolds-stress  tensor  by  hotwires  along  the  same 
streamlines.  However,  his  measurements  were  restricted  to  the  region  0.75  <  x 
<  0.95,  where  the  boundary  layer  was  thick  enough  for  probe  size  and  interfer¬ 
ence  effects  to  be  small. 
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(b)  Numerical  Grid 

As  in  the  previous  case,  the  solution  domain  is  defined  by  (0.3  <  x  < 
4.524,  rs  <  r  <  1.0),  and  covered  by  50  x  30  x  15  grid  points  in  the  axial, 
radial,  and  circumferential  directions,  respectively.  From  Section  IV. 2,  we 
recall  that  this  is  not  the  finest  gird  that  was  employed  to  study  the  grid 
dependence  of  the  present  method.  Some  views  of  the  numerical ly-generated 
grid  are  shown  in  Figure  20.  The  grid-control  functions  f1  used  to  obtain 
this  grid  were  given  in  equation  (41).  The  parameters  appearing  therein  are 
as  follows: 

0.10  if  A  <  0 

0.01  if  A  >  0 

d  =  0.2 

With  this  choice,  the  first  grid  node  just  off  the  hull  lies  in  the  range  40  < 
y+  <  250. 

(c)  Description  of  Results  and  Comparisons  with  Experiments 

We  have  already  mentioned  the  difficulties  of  interpolating  data  obtained 
in  boundary-layer  coordinates  to  make  comparisons  with  solutions  obtained  in 
numerical  coordinates  of  the  type  used  here.  We  shall  therefore  restrict  the 
comparisons  only  to  those  quantities  which  are  either  unaffected  by  the  choice 
of  coordinates  or  are  insensitive  to  differences  in  the  coordinates.  Among 
the  former  are  the  distributions  of  pressure  and  wall  shear  stress  on  the 
hull.  Figure  21  shows  the  distribution  of  pressure  and  friction  velocity  U. 
along  the  keel  and  the  waterline  planes  of  symmetry.  Girthwise  distributions 
of  pressure  and  the  magnitude  of  the  friction  coefficient  at  a  few  transverse 
sections  are  shown  in  Figure  22.  The  direction  of  the  wall  shear  stress  can 
be  seen  from  the  limiting  streamlines  plotted  in  Figure  23.  The  last  figure 
was  obtained  by  integration  of  the  wall  shear  stress  vectors  starting  from  the 
initial  section  at  x  =  0.3.  The  data  of  La^sson  (1974),  with  and  without 
blockage  correction  in  the  case  of  pressure,  are  also  shown  in  Figures  21  and 
22. 
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Figures  21  and  22  indicate  that  the  effect  of  wind  tunnel  blockage  is 
quite  significant  in  the  pressure  distributions  and  that  the  calculations  are 
in  quite  good  agreement  with  the  data  corrected  for  blockage.  It  is  seen  that 
at  midships  the  pressure  is  essentially  uniform  around  the  girth,  but  girth- 
wise  gradients  develop  as  the  geometry  changes.  These  gradients  become  large 
near  the  keel  as  the  transverse  surface  curvature  increases,  and  a  region  of 
low  pressure  develops  some  distance  from  the  keel.  The  associated  girthwise 
pressure  gradient  is  such  that  there  is  a  divergence  of  flow  out  of  the  keel 
and  a  thinning  of  the  boundary  layer  along  the  keel.  There  is  also  a  second¬ 
ary  motion  towards  the  lew  pressure  region  from  above.  This  convergence  of 
flow  from  both  sides  into  a  region  around  midgirth  leads,  as  we  shall  see,  to 
a  thickening  of  the  viscous  layer  there.  The  relationship  between  the  pres¬ 
sure  distribution  and  the  wall  shear  stress  is  not  immediately  obvious  from 
Figures  21  through  23,  but,  as  will  become  clear  later,  the  low  values  of  the 
wall  shear  stress  occur  where  the  viscous  layer  becomes  thick.  Figure  23 
shows  the  divergence  of  wall  streamlines  out  of  the  keel  over  much  of  the 
hull,  but  a  rather  complex  pattern  develops  over  the  stern.  The  strong  con¬ 
vergence  of  streamlines  in  an  area  just  above  the  keel  at  the  stern  is  fre¬ 
quently  interpreted  as  evidence  of  a  longitudinal  vortex  and  even  separ¬ 
ation.  However,  as  will  become  clear  as  we  examine  the  details  of  the  flow 
for  other  hulls,  this  pattern  by  itself  is  not  an  indication  for  either  a 
vortex  or  flow  separation. 

The  flow  pattern  on  the  waterplane  and  in  the  vertical  centerplane, 
constructed  by  interpolation  of  the  numerical  results,  is  shown  in  Figure 
24.  The  thickening  of  the  boundary  layer  over  the  stern  in  the  waterplane, 
and  the  thinning  of  the  boundary  layer  along  the  keel  are  clearly  seen  from 
these  views  as  is  the  evolution  of  the  three-dimensional  wake.  Contours  of 
the  axial  velocity  ( U)  and  velocity  vectors  projected  into  transverse  sections 
(V  and  W  components)  are  shewn  in  Figure  25  for  some  representative  sections 
in  the  range  0.5  <  x  <  1.1.  The  transverse  velocity  plots  also  give  an  indi¬ 
cation  of  the  grid  distribution.  At  midships,  the  secondary  flow  is  directed 
away  from  the  keel  as  well  as  from  the  waterplane.  This  leads  to  the  thicken¬ 
ing  of  the  layer  around  midgirth  which  is  already  evident  by  x  =  0.7.  This 
process  continues  as  the  secondary  motion  becomes  stronger,  particularly  near 
the  keel  where  the  hull  geometry  changes  more  rapidly.  We  now  see  that  the 
regions  of  low  wall  shear  stress  in  Figure  22  coincide  with  the  regions  of 
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thickening  of  the  viscous  layer.  By  x  =  0.95,  the  thickness  of  the  viscous 
layer  is  almost  as  large  as  the  draft  of  the  model.  At  and  downstream  of  the 
stern,  the  secondary  motion  gives  the  impression  of  a  longitudinal  vortex,  but 
the  strength  of  this  motion  diminishes  quite  rapidly. 

The  contours  of  the  axial  component  of  velocity  shown  in  Figure  25  may  be 
compared  with  those  of  the  magnitude  of  mean  velocity,  reconstructed  by 
Broberg  and  Larsson  (1984)  from  the  measurements  of  Larsson  (1974)  and  Lofdahl 
(1982),  which  are  shown  in  Figure  26.  We  note  again  not  only  the  difference 
in  the  quantities  but  also  the  fact  that  the  experimental  contours  are  not 
based  on  measurements  in  transverse  sections.  The  error  involved  is  small, 
however,  and  there  is  good  agreement  between  the  calculations  and  measurements 
with  regard  to  the  major  features  of  the  stern  flow. 

More  detailed  comparisons  between  experiments  and  calculations  are 
.attempted  in  Figure  27.  The  calculated  results  were  interpolated  to  determine 
the  variation,  with  distance  N  normal  to  the  transverse  section  of  the  hull, 
of  the  resultant  velocity  parallel  to  the  hull  (Q),  its  component  (C)  normal 
to  the  direction  of  the  calculated  streamline  at  the  location  of  the  measured 
boundary-layer  edge,  and  the  turbulent  kinetic  energy  (k),  at  each  of  the 
points  on  the  hull  for  which  experimental  data  are  available.  All  of  these 
quantities  are  relatively  insensitive  to  the  remaining  differences  between  the 
computational  and  experimental  locations.  Also,  we  note  that  C  is  the  so- 
called  crossflow  velocity  in  the  terminology  of  boundary-layer  theory. 

Figure  27a  shows  the  profiles  of  Q  and  C.  Figure  27b  shows  semi  logarith¬ 
mic  plots  of  Q  versus  N*  =  UQN/v,  the  format  of  the  Clauser  plots  which  were 
used  for  the  experimental  determination  of  the  wall  shear  stress  from  the 
measured  velocity  profiles.  Figure  27c  shows  the  distributions  of  the  turbu¬ 
lent  kinetic  energy.  Comparisons  are  presented  for  five  transverse  sections, 
namely,  x  =  0.50,  0.75,  0.85,  0.90,  and  0.95,  the  last  four  of  which  coincide 
with  those  where  the  turbulence  measurements  of  Lofdahl  were  made.  For  con¬ 
venience,  the  girthwise  locations  of  the  measurement  stations,  the  normals  to 
the  transverse  section  at  those  stations,  and  the  hull  section  just  downstream 
of  the  measurement  section  (from  Figure  19)  are  shown  in  each  figure.  The 
line  numbers  correspond  to  the  potential-flow  streamlines  calculated  by 
Larsson,  lines  1  and  10  being  the  keel  and  waterline,  respectively.  The  data 
are  shown  by  symbols  at  all  stations  where  measurements  were  made. 
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Considering  the  midship  section  first,  we  see  that  the  thickness  of  the 
boundary  layer  and  the  longitudinal  velocity  profiles  are  predicted  quite 
well.  The  crossflow,  which  is  small,  is  captured  only  qualitatively.  The 
Clauser  plots  indicate  that  the  calculations  reproduce  the  logarithmic  layer 
indicating  that  the  placement  of  the  near-wall  nodes  in  the  numerical  grid  is 
satisfactory  for  the  wall  functions  approach  to  be  applicable.  It  is  also 
seen  that  further  grid  refinement  can  be  made  by  pulling  the  near-wall  nodes 
closer  to  the  wall,  but  as  we  have  already  shown  in  Section  IV. 2,  adding 
points  in  the  grid  does  not  by  itself  provide  greater  accuracy.  More  impor¬ 
tantly,  the  agreement  between  the  calculations  and  data  in  the  inner  region  is 
indicative  of  the  accuracy  with  which  the  large  gradients  in  the  wall  layer 
are  resolved  in  the  numerical  solutions.  The  agreement  in  the  Clauser  plots 
is,  of  course,  reflected  in  the  agreement  shown  earlier  with  respect  to  the 
wall  shear  stress.  Finally,  there  is  nothing  remarkable  about  the  turbulent 
kinetic  energy  profiles  at  midships;  they  resemble  those  in  a  flat-plate 
boundary  layer. 

The  section  x  =  0.75  may  be  the  most  downstream  section  at  which  the  flow 
could  still  be  regarded  as  of  the  boundary-layer  type.  Here,  the  velocity 
magnitude  and  the  crossflow  component  are  predicted  with  reasonable  accuracy 
everywhere  except  at  line  6,  where  it  is  the  data  which  appear  to  be  spur¬ 
ious.  The  profiles  of  turbulent  kinetic  energy  at  this  section  begin  to  show 
a  trend  which  becomes  very  pronounced  further  downstream.  We  observe  that  the 
calculated  values  are  generally  larger  than  the  measurements  and,  at  line  5, 
which  is  roughly  where  the  thickness  of  the  viscous  layer  is  the  greatest,  the 
k  distribution  develops  a  plateau. 

The  final  three  sections  considered  in  Figure  27  (x  =  0.85,  0.90,  0.95) 
lie  in  the  stern-flow  region  where  the  viscous  layer  can  no  longer  be  regarded 
as  a  boundary  layer.  From  the  results  at  these  sections  we  make  the  following 
observations.  First,  the  variation  of  the  thickness  of  the  viscous  layer, 
around  the  girth  and  in  the  axial  direction,  are  predicted  with  considerable 
accuracy.  Secondly,  the  mean  velocity  profiles  are  described  quite  well  by 
the  calculations  although,  as  evidenced  by  the  logarithmic  plots,  there  is  a 
tendency  for  the  calculated  velocities  to  be  somewhat  larger  in  the  near  wall 
region.  Thirdly,  the  general  shapes  of  the  crossflow  profiles  are  predicted 
but  there  exist  differences  in  magnitude.  The  differences  may,  however,  be  of 
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the  same  order  as  the  uncertainties  of  interpolation.  Fourthly,  the  predicted 
turbulent  kinetic  energy  is  considerably  larger  than  that  measured.  Finally, 
we  note  a  very  characteristic  two-layer  feature  developing  around  midgirth 
(see  line  5  at  x  =  0.85;  and  lines  3,5,7  at  x  =  0.90  and  0.95)  which  suggests 
the  existence  of  a  thin  layer  of  fluid  close  to  the  wall  in  which  k  diminishes 
rapidly,  and  a  much  larger  layer  farther  out  where  there  is  a  plateau  in  k  and 
a  gradual  decrease  to  zero  outside  the  viscous  region.  Although  there  is  a 
hint  of  such  a  two-layer  structure  in  the  calculated  profiles,  it  is  obvious 
that  the  present  turbulence  model  fails  to  capture  the  changes  that  are  taking 
place  in  the  turbulence  in  the  flow  over  the  stern. 

(d)  Summary 

The  calculations  for  the  SSPA  liner  confirm  the  capability  of  the  present 
numerical  method  to  predict  the  essential  features  of  the  pressure  and  mean- 
velocity  fields  in  considerable  detail.  This  is  particularly  surprising  in 
view  of  the  differences  between  the  calculated  and  measured  turbulent  kinetic- 
energy  distributions,  and  leads  to  the  conclusion  that  the  turbulence  in  the 
outer  part  of  the  viscous  layer  is  inactive  and  is  simply  being  convected  with 
the  mean  flow.  The  success  of  the  calculations  must  therefore  come  from  the 
use  of  the  wall  functions  which  essentially  reproduce  the  observed  logarithmic 
layer  close  to  the  surface.  This  also  suggests  that  further  improvements  are 
most  likely  to  result  from  refinement  of  the  treatment  of  the  flow  in  the 
near-wall  region.  Since  the  experiments  on  this  hull  did  not  extend  into  the 
wake,  it  is  not  possible  to  comment  on  the  quality  of  the  wake  calculations. 

Finally,  we  note  that  these  calculations  have  also  pointed  out  the  impor¬ 
tance  of  wind-tunnel  blockage,  particularly  with  regard  to  the  hull  pressure 
distribution.  The  effects  of  blockage  on  the  velocity  distributions  are  too 
small  to  be  seen  on  the  various  plots.  Because  the  present  method  can  also 
take  into  account  wind-tunnel  walls  by  appropriate  relocation  of  the  external 
boundary  of  the  solution  domain,  it  would  be  of  interest  to  carry  out  such  a 
calculation  to  study  the  blockage  effect  in  greater  detail. 


V.3  HSVA  TANKER 


(a)  Hull  Geometry  and  Available  Experimental  Information 

The  offsets  of  the  HSVA  Tanker  are  shown  in  Figure  28.  The  boundary 
layer  measurements  by  Hoffman  (1976),  and  the  detailed  stern-flow  measurements 
of  Wieghardt  and  Kux  (1980)  and  Wieghardt  (1982,  1983),  were  made  on  a  double 
model  of  this  hull  in  a  1.2  m  square,  slotted-wall  wind  tunnel,  in  which  the 
turbulence  level  was  of  the  order  of  1  percent.  The  slotted  walls  are  usually 
meant  to  reduce  or  eliminate  blockage  effects.  The  model  was  supported  in  the 
tunnel  by  means  of  wires  and  a  sting  at  the  stern.  The  nominal  length  of  the 
model  was  2.74  m  but  for  reference  length  we  have  used  the  length  between 
perpendicuars,  L  =  2.6G4  m.  The  model  beam  and  draft  are  0.43  m  and  0.15  m, 
respectively,  and  the  block  coefficient  is  0.85. 

The  locations  of  the  measurement  sections  in  the  experiments  of  Wieghardt 
and  Kux  are  shown  in  Figure  29.  The  data  were  obtained  at  nine  axial  sections 
in  the  range  -157  mm  <  Y  <  200  mm,  in  the  notation  of  Figure  29.  This  corre¬ 
sponds  to  0.9103  <  x  <  1.044  in  the  present  notation.  The  measurements  were 
made  by  means  of  a  five-hole  pitot  probe,  and  included  the  three  components  of 
mean  velocity  (U,V,W)  in  Cartesian  coordinates  and  pressure.  The  closely 

spaced  measurements  enabled  the  determination  of  the  mean  vorticity  field. 

_  f) 

The  calculations  were  performed  for  a  Reynolds  number  of  5  x  10  which  corre¬ 
sponds  to  the  experimental  value. 

(b)  Numerical  Grid 

As  in  the  two  previous  cases,  the  chosen  solution  domain  U  defined  by 
0.3  <  x  <  4.524,  r$  <  r  <  1.0  and  was  covered  by  50  x  30  x  15  grid  points  in 
the  axial,  radial,  and  circumferential  airecticns,  respectively.  Some  views 
of  the  numerical ly-generated  gr^d  are  shown  in  Figure  30.  The  grid-control 
functions  f1  used  to  obtain  this  grid  were  given  in  equation  (41).  The  param¬ 
eters  appearing  therein  are  as  follows: 

0.20  if  A  <  0 

0.02  if  A  >  0 

d  =  0.25 
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With  this  choice,  the  first  grid  node  just  off  the  hull  lies  in  the  range  50  < 
y+  <  250. 

(c)  Description  of  Results  and  Comparisons  with  Experiments 

Figure  31  shows  the  distribution  of  pressure  and  friction  velocity  along 
the  keel  and  the  waterline  planes  of  symmetry,  and  the  girthwise  distributions 
of  pressure  and  the  magnitude  of  the  friction  coefficient  are  shown  in  Figure 
32  at  a  few  transverse  sections.  In  these  figures  we  have  included  the  data 
from  the  boundary- layer  experiments  of  Hoffman  (1976)  at  a  Reynolds  number  of 
4.8  x  10®.  The  section  x  =  0.9418  corresponds  to  the  most  downstream  measure¬ 
ment  section  in  his  experiments  and  to  the  J  =  -73  mm  section  in  the  later 
experiments  of  Wieghardt  and  Kux.  From  Figure  31  we  see  that  the  pressure 
distribution  is  in  good  agreement  with  the  data  except  along  the  keel  at  the 
stern.  The  lower  predicted  values  of  friction  velocity  along  the  keel  are 
consistent  with  the  differences  in  the  pressure  distributions.  These  differ¬ 
ences  are  again  evident  from  the  girthwise  variations  at  x  =  0.9417  shown  in 
Figure  32.  These  suggest  that  the  stern  flow  in  the  neighborhood  of  the  keel 
is  not  being  described  with  sufficient  accuracy  This  may  be  due,  at  least  in 
part,  to  differences  between  the  actual  hull  geometry  and  that  implied  in  the 
calculations  by  the  chosen  numerical  grid. 

The  direction  of  the  wall  shear  stress  can  be  seen  from  the  limiting 
streamlines  plotted  in  Figure  33.  This  shows  a  region  around  midgirth  at  the 
stern  into  which  the  wall  streamlines  converge.  This  is  the  area  of  maximum 
thickness  of  the  viscous  layer.  These  features,  like  those  of  the  results 
presented  in  the  previous  two  figures,  are  qualitatively  similar  to  those 
discussed  earlier  for  the  SSPA  liner. 

The  rather  complex  stern  geometry  of  the  tanker  and  some  general  features 
of  the  resulting  flow  pattern  are  seen  from  Figure  34.  Comparison  with  the 
corresponding  figure  (Figure  24)  for  the  more  slender  SSPA  hull  shows  that  the 
flow  pattern  is  qualitatively  similar.  However,  these  views  of  the  flow  tend 
to  mask  the  very  significant  differences  that  exist. 

Figure  35  presents  comparisons  between  the  experimental  data  and  calcula¬ 
tions  with  respect  to  the  mean-velocity  field.  The  contours  of  axial  velocity 
and  projections  of  the  velocity  vectors  in  transverse  sections  are  shown  at 
six  sections.  Particularly  noteworthy  here  is  the  very  fine  grid  used  for  the 


measurements.  Here,  the  measurements  and  calculations  correspond  to  trans¬ 
verse  sections  and  therefore  the  comparisons  do  not  involve  the  uncertainties 
which  were  present  in  the  case  of  the  SSPA  liner. 

At  the  most  upstream  section,  x  =  0.6433,  we  see  that  the  boundary  layer 
is  thin  except  around  the  turn  of  the  bilge.  Although  the  local  crossflow  is 
quite  weak,  the  thicker  boundary  layer  is  the  result  of  flow  convergence 
towards  the  turn  over  some  distance  upstream.  By  the  next  section,  x  = 
0.9103,  the  boundary  layer  has  thickened  considerably  near  midgirth.  The 
following  section,  x  =  0.9418,  corresponds  to  the  last  measurement  section  in 
the  experiments  of  Hoffman  and  the  data  shown  in  Figure  lb.  Comparison  of  the 
calculations  with  the  recent  and  more  detailed  data  indicates  that  the  general 
features  of  the  U  contours  and  the  transverse  velocity  components  are  pre¬ 
dicted  at  all  sections  but  there  are  differences  particularly  close  to  the 
hull  around  z  =  -0.04.  The  gaps  in  the  data  at  these  depths  at  x  =  0.9549  and 
0.9692  are  regions  of  low  velocity  which  could  not  be  resolved  by  the  pitot 
probe.  At  the  three  downstream  sections,  some  of  the  velocity  contours  are 
missing  partly  because  the  predicted  velocities  are  larger  and  partly  due  to 
the  use  of  the  wall  functions.  The  wall  functions  could  of  course,  be  used  to 
reconstruct  the  contours  close  the  hull  and  wake  centerplane  but  this  has  not 
been  done  here.  Comparisons  of  specific  contours,  for  example  U  =  0.5  and  0.9 
at  the  last  three  sections  reveal  that  the  predictions  are  in  remarkably  good 
agreement  with  the  measurements.  The  solutions  also  capture  the  major  fea¬ 
tures  of  the  secondary  motion,  including  the  magnitudes,  but  this  is  not 
immediately  obvious  from  Figure  35. 

In  order  to  make  a  detailed  assessment  of  the  predictions,  the  results  of 
the  calculations  were  interpolated  to  obtain  the  distributions  of  the  three 
velocity  components  ( U , V , W )  and  pressure  (p)  in  the  Cartesian  coordinates 
which  were  employed  for  the  measurements.  This  makes  it  possible  to  plot  the 
profiles  of  these  parameters  in,  say,  the  horizontal  direction  (y)  at  differ¬ 
ent  depths  (z  =  constant).  The  results  are  shown  in  Figure  36.  Although 
there  is  a  great  deal  of  similarity  in  the  results  because  of  the  closely 
spaced  stations,  we  have  included  all  ten  transverse  sections  here  to  facil¬ 
itate  future  reference  and  compar'sons  by  others  using  different  methods.  In 
general,  we  see  that  the  calculations  are  in  rather  remarkable  agreement  with 
the  data  in  all  respects  except  in  small  regions  close  to  the  hull  or  in  the 
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neighborhood  of  the  wake  centerplane  just  downstream  of  the  ship.  Before 
discussing  these  features  it  is  useful  to  point  out  that  y  in  these  profiles 
is  measured  from  the  ship  centerplane  and  therefore  one  gets  a  rough  idea  of 
the  hull  shape  at  ea^h  section.  Also,  z  =  0  is  the  waterplane  and  the  keel  is 
at  z  =  -0.051. 

Reference  to  Figure  29  shows  that  the  first  five  sections  are  located  in 
the  region  where  the  keel  is  still  horizontal  (x  =  0.9621  corresponding  to 
section  X  =  -19  mm),  the  next  section  is  just  where  the  stern  is  almost  per¬ 
pendicular  (x  =  0.9692  corresponds  to  X  =  0),  and  the  remaining  four  sections 
are  located  downstream  of  the  propeller  plane.  It  is  important  to  keep  in 
mind  the  local  geometry,  and  particularly  the  changes  in  it,  in  order  to 
understand  the  results  shown  in  Figure  36.  It  is  convenient  to  discuss  these 
in  three  stages,  roughly  corresponding  to  the  three  regions  identified  above. 

In  the  first  region,  upto  x  =  0.9621  (Figure  36a-e),  we  see  that  the  U 
and  V  components  are  in  good  agreement  with  the  data  but  the  vertical  compon¬ 
ent  W  shows  a  characteristic  difference  that  persists  downstream.  The  magni¬ 
tude  of  the  W  component  changes  rather  rapidly  with  distance  from  the  hull 
surface  in  the  region  where  the  surface  is  nearly  vertical,  i.e.,  where  W  is 
nearly  parallel  to  the  surface.  At  the  two  downstream  sections,  we  also 
observe  increasing  differences  in  the  horizontal  V  component.  The  prediction 
of  the  pressure  field  is  quite  satisfactory,  except  in  small  layers  close  to 
the  hull  which  roughly  coincide  with  the  regions  in  which  we  see  differences 
in  the  W  velocity  component.  There  are  two  possible  explanations  for  the 
differences  between  the  calculations  and  measurements.  One  is  that  the  wall- 
functions  approach  employed  here  coir  not  completely  capture  the  details  of 
the  secondary  motion.  The  second  a  "h  more  likely  possiblity  is  that  the 
numerical  grid  does  not  resolve  the  local  geometry  in  sufficiently  fine  detail 
to  predict  the  secondary  motions  and  the  associated  pressure  field. 

The  plots  for  x  =  0.9692  (Figure  36f)  rev*  1  that  the  differences  between 
the  calculations  and  measurements  now  extend  to  the  U  and  V  components,  and 
the  characteristic  features  of  the  W  component  noted  above  continue.  A  care¬ 
ful  examination  of  these  differences  indicates  that  they  stem  principally  from 
the  truncation  of  the  local  geometry  implied  by  the  numerical  grid  employed, 
and  from  the  way  in  which  the  hull  and  wake  boundary  conditions  are  enforced 
in  the  local  numerical  cells  at  this  section  within  the  framework  of  the 
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staggered  grid  used  in  the  present  velocity-pressure  coupling  method.  How¬ 
ever,  previous  experience  with  problems  involving  discontinuities  of  boundary 
conditions,  from  no  slip  on  the  hull  to  zero  stress  on  the  wake  centerplane, 
suggests  that  improved  results  can  be  obtained  by  more  refined  grids  and  a 
more  direct  treatment  of  the  flow  in  the  near-wall  layer  than  made  in  the 
wall-functions  approach. 

The  results  at  the  final  four  sections  (x  >  0.9786,  Figure  36g-j)  indi¬ 
cate  the  development  of  the  near  wake.  Here,  the  differences  which  were 
generated  upstream  continue.  Also,  we  see  that  the  calculations  predict  a 
rather  rapid  increase  in  axial  velocity  in  the  wake  centerplane,  a  feature 
which  was  also  observed  in  the  case  of  the  Wigley  hull  (see  Figure  18f-l). 
This  is  most  likely  due  to  the  use  of  the  wall-functions  which  do  not  resolve 
the  flow  in  the  near-wall  layer  on  the  hull  upstream.  It  is  interesting  to 
note  the  decay  of  the  calculated  as  well  as  the  measured  secondary  motion  with 
downstream  distance.  At  the  last  section  the  calculated  U  and  V  components 
are  again  in  reasonably  good  agreement  with  data  while  the  differences  in  W 
and  p  persist. 

As  we  have  already  noted,  the  measurements  of  Wieghardt  and  Kux  were  made 

in  such  detail  that  it  was  possible  to  carry  out  the  differentiations  required 

to  determine  the  mean  vorticity  field.  In  spite  of  the  uncertainties  involved 

in  processing  both  the  measured  and  computed  data,  a  comparison  between  the 

computed  and  measured  vorticity  is  a  particularly  acid  test  of  a  computational 

method  because  it  is  an  examination  of  differences  in  velocity  gradients 

rather  than  just  velocity  components.  Such  a  comparison  is  presented  in 

Figures  37  and  38.  Figure  37  shows  the  contours  of  the  longitudinal  component 

of  vo’  <.icity,  [3 W/3y  -  3V/3z],  while  Figure  38  shows  the  contours  of  the 

2  2  2  M2 

magnitude  of  the  vorticity  vector,  (<•>  +  &>  +  w  )  .  These  plots  show  that 

x  y  z 

the  calculations  reproduce  not  only  the  general  features  of  the  contours  but 
also  the  magnitudes.  The  differences  close  to  the  hull  and  around  the  wake 
centerplane  are  to  be  expected  from  the  differences  in  the  velocity 
distributions  discussed  above.  An  additional  and  new  feature  that  becomes 
apparent  from  the  vorticity  comparisons,  however,  is  that  the  calculations 
indicate  a  somewhat  larger  diffusivity  compared  with  the  measurements.  This 
is  presumably  due  to  an  overestimation  of  the  eddy  viscosity  by  the 
standard  k-e  turbulence  model.  We  have  already  commented  on  this  feature  in 


connection  with  the  prediction  of  the  turbulent  kinetic  energy  in  the  previous 


cases. 


(d)  Discussion 

The  HSVA  tanker  represents  a  more  severe  test  of  the  calculation  method 
than  the  previous  two  hulls  because  of  its  high  block  coefficient  and  rather 
abrupt  changes  in  stern  geometry.  In  addition  to  the  observations  that  have 
been  made  on  the  basis  of  the  previous  calculations,  the  comparisons  between 
the  nunerical  results  and  the  very  extensive  experimental  data  for  the  present 
case  suggest  the  need  for  (a)  more  refined  treatment  of  the  stern  geometry, 
(b)  better  resolution  of  the  flow  in  the  near-wall  layer,  and  (c)  some  correc¬ 
tion  for  the  higher  rates  of  diffusion  of  mean  vorticity  predicted  by  the 
calculations. 

V.4  SRI 07  ORE  CARRIER 

(a)  Hull  Geometry  and  Experimental  Information 

The  offsets  of  this  hull  are  shown  in  Figure  39.  Experiments  on  a  double 
model  of  this  form  were  conducted  in  the  1.8  m  x  1.8  m  wind  tunnel  of  Osaka 
University  (Okajima  et  al.,  1985  and  Nishio  et  al.,  1988).  The  2  m  long  model 
was  supported  in  the  tunnel  by  means  of  wires  along  the  keel.  The  model  beam 
and  draft  are  0.2927  m  and  0.1126  m,  respectively,  and  the  block  coefficient 
is  0.826.  The  data  were  obtained  at  five  axial  sections  in  the  range  0.8  <  x 
<  1.0  in  the  present  notation.  The  measurements  were  by  means  of  a  five-hole 
pitot  probe,  and  included  the  three  components  of  mean  velocity  (U,V,W)  in 
Cartesian  coordinates  and  pressure.  The  closely  spaced  measurements  enabled 
the  determination  of  the  mean  vorticity  field.  The  calculations  were  per¬ 
formed  for  a  Reynolds  number  of  2.4xl06  which  corresponds  to  the  experimental 
value. 

(b)  Numerical  Grid 

The  numerical  solution  domain  chosen  for  this  case  is  the  same  as  in  the 
previous  cases  and  is  defined  by  0.3  <  x  <  4.524,  r$  <  r  <  1.0  .  This  was 
again  covered  by  50  x  30  x  15  grid  points  in  the  axial,  radial,  and  circumfer¬ 
ential  directions,  respectively.  Some  views  of  the  numerically-generated  grid 
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are  shown  in  Figure  40.  The  grid-control  functions  f1  used  to  obtain  this 


grid  were 

given 

in  equation  (41).  The  parameters  appearing  therein  are  as 

follows: 

(0.10 

if  A  <  0 

c  = 

<0.01 

if  A  >  0 

d  = 

0.2 

With  this  choice,  the  first  grid  node  just  off  the  hull  lies  in  the  range  60  < 
y+  <  250. 

(c)  Description  of  Results  and  Comparisons  with  Experiments 

Comparison  of  Figure  39  with  Figure  28  shows  the  difference  in  the  bow 
and  stern  shapes  of  the  SR107  and  the  HSVA  hulls  although  both  are  relatively 
full  forms.  The  measurements  in  the  two  experiments  were  made  in  a  very 
similar  manner  and  therefore  it  is  convenient  to  present  the  results  for  the 
present  case  in  essentially  the  same  format.  Thus,  Figures  41  and  42  show  the 
hull  pressure  and  friction  distributions.  Figure  43  shows  the  calculated  wall 
streamlines,  Figure  44  shows  the  flow  pattern  in  the  waterplane  and  in  the 
vertical  centerplane,  and  Figures  45  through  48  give  details  of  the  velocity 
and  vorticity  fields.  In  general,  these  results  are  qualitatively  similar  to 
the  corresponding  results  for  the  HSVA  hull.  However,  the  level  of  agreement 
between  the  measured  and  computed  pressure  distributions  shown  in  Figures  41 
and  42  is  particularly  noteworthy. 

An  overview  of  the  velocity  rield  is  given  by  Figure  45  which  shows  the 
contours  of  the  axial  velocity  component  and  projections  of  the  velocity 
vectors  in  five  transverse  sections.  These  may  be  compared  with  Figure  35. 
The  wake  of  the  support  wires  at  the  keel  is  seen  from  the  results  at  x  =  0.8, 
and  this  is  responsible  for  the  somewhat  thicker  boundary  layer  measured  along 
the  keel  further  downstream.  The  reason  for  the  abrupt  thickening  of  the 
boundary  layer  around  the  turn  of  the  bilge  shown  by  the  data  at  x  =  0.8  is 
not  clear  but  its  effect  is  observed  at  almost  a'l  stations  downstream.  The 
difference  between  the  calculations  and  experiments  in  this  regard  is  most 
likely  due  to  the  inadequacy  of  the  simple  initial  conditions  used  in  the 
calculations  at  x  =  0.3.  Prestmably,  the  bulbous  bow  of  this  model  induces  a 
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level  of  three  dimensionality  that  is:  not  captured  by  the  initial  condi¬ 
tions.  The  initial  conditions  ccdld  of' course*  have- teen  adjusted' to  obtain  a 
better  fit  to  the  measurements  at  the- "First  measurement  section  but  this  has 
not  been  attempted  here. 

The  results  of  Figure  45  show  yet  another  feature  that  was  absent  in  the 
previous  cases.  The  SR107  model  included'  a  propeller  boss  and  a  rudder  post, 
which  are  visible  in  the  plots  for  x  =  0.975.  These  details,  however,  were 
not  considered  in  the  present  calculations  and  the  stern  sections  were  simply 
faired  as  shown.  It  is  clear  from  the  comparisons  between  the  calculations 
and  measurements,  particularly  at  sections  x  =  0.95,  0.975,  and  1.00,  that 
there  is  a  considerable  local  effect  of  these  features  on  the  velocity 
field.  The  differences  are  of  course  more  pronounced  in  the  transverse  com¬ 
ponents. 

The  aforementioned  differences  should  be  kept  in  mind  in  assessing  the 
detailed  comparisons  of  the  velocity  components  and  pressure  presented  in 
Figure  46.  The  format  is  the  same  as  that  used  in  the  case  of  the  HSVA  tanker 
(see  Figure  36).  It  is  quite  remarkable  that  the  predicted  values  of  all 
three  velocity  components  are  in  good  agreement  with  the  data  everywhere  at 
the  first  section,  x  =  0.8,  except  around  z  =  -0.045,  a  region  that  corre¬ 
sponds  to  the  thick  boundary  layer  mentioned  earlier.  The  pressures  measured 
at  this  section  show  a  rather  curious  behavior  although  the  calculations  and 
data  are  in  good  agreement  very  close  to  the  hull.  In  fact,  the  measured 
trends  are  rather  difficult  to  understand.  The  persistence  of  this  feature  at 
other  stations  downstream  would  suggest  a  consistent  error  in  the  data. 

The  distributions  of  the  velocity  components  at  x  =  0.9  are  very  similar 
to  those  at  x  =  0.8.  Some  differences  begin  to  appear  at  x  =  0.95.  These  are 
consistent  with  the  effects  of  the  propeller  boss  and  rudder  post  mentioned 
above.  A  particularly  noteworthy  feature  of  the  results  at  x  =  0.95  and  1.00, 
however,  is  their  marked  resemblance  with  those  obtained  around  the  stern  of 
the  HSVA  tanker.  Comparison  of  the  SR107  results  at  x  =  1.00  with  those  for 

i 

the  HSVA  tanker  at  x  =  0.988  (Figure  36),  for  example,  reveals  very  similar 
distributions  of  all  velocity  components.  This  leads  us  to  reiterate  the 
previous  observations  concerning  the  need  for  a  better  resolution  of  the  local 
geometry  and  flow. 


Contours  of  the  longitudinal  component  of  vorticity,  o>  =  (3 W/3y-3 V/3 z) , 

x  2 

are  shown  in  Figure  47,  and  Figure  48  shows  contours  of  {(3(J/3z) 

o  i/p  221/2 

+  (3U/3y)  )  ,  which  is  approximately  )  ,  the  resultant  of  the 

y-and  z-components  of  vorticity.  It  is  clear  that  the  calculations  are  in 
remarkable  general  agreement  with  the  experimental  results.  Some  of  the 
details  are  not  completely  captured  by  the  solutions  for  reasons  which  we  have 
already  discussed.  It  is,  however,  interesting  to  note  that  the  overall 
vorticity  field  is  not  markedly  affected  by  the  presence  of  the  propeller  boss 
and  the  rudder  post.  This  is  evidently  due  to  the  fact  that  while  they 
produce  significant  distortion  of  the  velocity  and  pressure  fields,  they  do 
not  generate  large  amounts  of  vorticity  that  can  be  seen  outside  their 
immediate  wakes. 

(d)  Oiscussion 

The  solutions  for  the  SR107  hull  and  comparisons  with  the  corresponding 
experiments  have  brought  forth  three  issues  which  were  not  evident  from  the 
previous  test  cases.  First,  we  observe  that  a  hull  with  a  bulbous  bow  may 
require  a  more  careful  treatment  of  the  initial  conditions  than  has  been 
employed  in  the  present  calculations.  However,  the  differences  observed  in 
this  respect  may  have  been  exaggerated  in  the  present  case  by  the  rather  low 
Reynolds  number  of  the  experiments  and  perhaps  an  ineffective  transition 
device.  Secondly,  the  local  distortions  of  the  velocity  and  pressure  field 
introduced  by  such  practical  geometrical  features  as  a  propeller  boss  and  a 
rudder  post  need  to  be  taken  into  account  in  comparisons  between  experiments 
and  calculations.  Incorporation  of  such  geometrical  features  would  require 
further  development  of  certain  aspects  of  the  computational  method.  Thirdly, 
the  very  marked  similarity  of  the  results  for  the  SRI 07  and  HSVA  hulls  points 
to  the  need  for  refinement  of  the  treatment  of  the  near-wal1  region  in  stern 
flows. 

V.5  SERIES  60,  <}>  =  0.60 

(a)  Hull  Geometry  and  Experimental  Information 

The  geometry  of  this  hull  is  well  known  and  is  shown  in  Figure  49.  The 
model  used  in  the  experiments  was  fitted  with  a  propeller  boss  (stern  tube)  as 


shown  in  Figure  50  but,  as  in  the  previous  case,  the  calculations  correspond 
to  the  bare  hull  shown  in  Figure  49.  Also,  as  noted  earlier,  these 

experiments  differ  from  the  others  in  that  they  were  conducted  in  a  towing 

tank.  A  detailed  description  of  the  experiments  is-  given  by  Toda  et  al., 

(1988).  For  the  purposes  of  the  present  work,  it  suffices  to  note  the 

following.  The  measurements  were  made  on  4  m  long  models  at  a  Froude  number 
Fn  =  0.16  and  a  Reynolds  nimber  Re  =  3.2  x  10®.  Two  different  models  were 
employed,  one  for  the  hull  pressure  distributions  and  the  Other  for  the  flow- 
field  measurements.  The  calculations  have  been  made  with  a  flat  free  surface 
(Fn  =  0)  and  plane-of-symmetry  boundary  conditions  on  it.  The  available  data 
include  the  hull  pressure  distribution  measured  by  surface  pressure  taps,  and 
profiles  of  the  velocity  components  (U,V,W)  and  pressure  p  measured  by  a  five- 
hole  pitot  probe.  The  latter  data  were  obtained  at  numerous  transverse 
sections  in  the  range  0.50  <  x  <  1.10.  The  use  of  Cartesian  coordinates,  as 
in  the  experiments  on  the  HSVA  and  the  SR107  hulls,  greatly  facilitates 
detailed  comparisons  between  calculations  and  experiments. 

(b)  Numerical  Grid 

The  numerical  solution  domain  chosen  for  this  case  is  the  same  as  in  the 
previous  cases  and  is  defined  by  0.3  <  x  <  4.524,  r$  <  r  <  1.0  .  This  was 
again  covered  by  50  x  30  x  15  grid  points  in  the  axial,  radial,  and  circumfer¬ 
ential  directions,  respectively.  Some  views  of  the  numerical ly-generated  grid 
are  shown  in  Figure  51.  The  grid-control  functions  f1  used  to  obtain  this 
grid  were  given  in  equation  (41).  The  parameters  appearing  therein  are  as 
follows: 

(0.10  if  A  <  0 
c  =  ] 

(0.01  if  A  >  0 
d  =  0.2 

With  this  choice,  the  first  grid  node  just  off  the  hull  lies  in  the  range  50  < 
y+  <  250. 
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|  (c)  Description  of  Results  and  Comparisons  with  Experiments 

j 

j  The  variation  of  pressure  and  friction  velocity  along  the  keel  and  the 
j  waterline  are  shown  in  Figure  52.  No  data  is  shown  for  the  friction  velocity 

because  the  velocity  profiles  have  not  yet  b«:en  analyzed  through  Clauser 
plots.  On  the.  other  hand,-  two  sets  of  data  are  shown  for  the  pressure.  These 
were  obtained  during  the  pitot  measurements  on  the  port  and  starboard  sides. 
The  data  on  the  waterline  are  those  measured  just  below  the  free  surface.  The 
close  agreement  of  the  calculations  with  the  measured  pressure  distribution 
along  the  waterline  indicates  that  the  effect  of  the  free  surface  at  the  low 
Froude  number  of  these  tests  is  negligible. 

>  The  girthwise  distributions  of  pressure,  measured  by  means  of  surface 

pressure  taps,  at  four  typical  transverse  sections  are  shown  in  Figure  53. 
Also  shown  there,  for  future  reference,  are  the  distributions  of  the  friction 
coefficient  Cf.  The  major  features  of  the  hull  pressure  distribution  are 
predicted  quite  well  except  that  there  appears  to  be  a  nearly  constant  differ¬ 
ence  between  the  calculations  and  data  at  the  first  three  sections,  and  a 
disagreement  in  a  region  above  the  keel  at  x  =  0.975.  The  former  could  be  due 
to  a  difference  between  the  reference  pressure  used  in  the  experiments  and  the 
true  freestream  pressure.  The  latter,  on  the  other  hand,  is  associated  with 
the  presence  of  the  propeller  boss  in  the  experiments. 

Comparison  of  the  results  for  the  Series  60  with  those  obtained  for  the 
SSPA  liner  indicates  that  the  general  features  of  the  pressure  and  friction 
distributions  on  these  hulls  are  quite  similar.  This  is  also  the  case  for  the 
wall  streamlines  shown  in  Figure  54,  and  the  flow  patterns  in  the  waterplane 
and  the  vertical  centerplane  shown  in  Figure  55. 

The  contours  of  constant  axial  velocity  and  the  transverse  components  of 
velocity  are  shown  in  Figure  56.  It  is  clear  from  the  plots  of  the  data  that 
the  measurements  could  not  be  made  very  close  to  the  free  surface  which  was 
relatively  flat  over  the  stern  region  and  coincided  with  the  load  waterline,  z 
=  0.  Excellent  agreement  is  observed  between  the  calculated  and  measured  ax¬ 
ial  velocity  contours  upto  about  x  =  0.95.  A  close  examination  of  the  trans¬ 
verse  velocity  field  indicates  that  some  differences  begin  to  appear  above  and 
below  z  =  -0.032  which  is  where  the  centerline  of  the  propeller  boss  is  locat¬ 
ed.  This  is  undoubtedly  an  upstream  influence  of  the  boss.  The  results  for  x 
=  0.975  and  downstream  of  there  clearly  show  the  boss  and  its  effect.  It  is 
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seen  that  this  influence  is  rather  localized,  as  was  the  case  for  the  SR107 
model,  and  the  agreement  between  the  predictions  and  experiments  continues  to 
be  quite  good  elsewhere. 

Finally,  detailed  distributions  of  the  velocity  components  (U,V,W)  and 
pressure  p,  are  shown  in  Figure  57  in  a  format  that  should  be  quite  familiar 
by  now.  In  view  of  the  discussion  that  has  already  been  made  of  the  results 
for  the  HSVA  tanker  (Figure  36)  and  the  SR107  Ore  Carrier  (Figure  46),  we 
shall  point  out  only  the  two  major  distinguishing  features  of  these  results. 
First  of  all,  we  see  a  high  level  of  agreement  between  the  calculations  and 
experimental  data  in  almost  all  respects  upto  x  =  0.9.  There  are  occasional 
differences  which  appear  to  be  due  to  scatter  in  the  data  rather  than  to  some 
flow  feature  that  is  not  being  resolved  by  the  calculations.  Secondly,  the 
regions  of  disagreements  at  the  downstream  stations  clearly  define  the  influ¬ 
ence  of  the  local  change  in  the  stern  geometry  made  by  the  propeller  boss. 
The  fact  that  this  influence  is  localized  in  the  transverse  sections  is  seen 
from  the  continued  agreement  between  the  experimental  data  and  the  bare-hull 
calculations. 

(d)  Discussion 

These  calculations  have  revealed  the  similarities  in  the  stern  flow  of 
the  Series  60  and  the  SSPA  liner.  They  have  also  indicated  that  in  the 
absence  of  a  wave  system  the  effect  of  the  free  surface  on  the  viscous  flow 
appears  to  be  quite  weak.  The  detailed  *  wparisons  in  the  velocity  profiles 
clearly  delineate  the  effects  of  local  hull  modifications  such  as  that  of  a 
propeller  boss.  These  comparisons  are  presented  here  in  detail  because  they 
will  be  useful  in  the  assessment  of  future  calculations  that  include  such 
complexities. 


VI.  SUMMARY  AND  CONCLUSIONS 

A  f ul ly-el 1 iptic  numerical  method  for  the  solution  of  the  complete 
Reynold-averaged  Navier  Stokes  equations  for  general  three-dimensional  flows 
has  been  presented.  The  method  uses  numerically-generated  nonorthogonal 
coordinates  while  retaining  the  velocity  components  in  a  convenient  orthogonal 
system.  For  turbulent  flows,  closure  of  the  equations  is  effected  by  the  well 
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known  k-e  turbulence  model  with  the  generally  accepted  values  of  the  constants 
in  the  model'  and  a  two-point  wall-function  approach.  The  various  numerical 
features- of  the  method  have  been  described  and  evaluated  in  detail. 

The  numerical  method  described  here  is  built  into  a  computer  code 
RANSTERN  which  consists  of  three  programs.  The  first  of  these  maps  the 
physical  hull  surface  into-  the  rectangular  plane-n  =  1  in  the  transformed 
domain  (see  Figures  2  and  3).  Spline  functions  are  used  for  interpolation  of 
the  numerical  grid  on  the  hull  surface.  The  surface  grid  thus  obtained  is 
used  to  provide  the  boundary  conditions  for  the  second  program  which  generates 
the  body-fitted  coordinate  system  for  the  whole  computational  domain  as 
described  in  Section  III. 2. (a),  the  flow  equations  are  solved  in  the  third 
program  using  the  methods  described  in  Section  1 1 1 .2. (b) .  For  the 
calculations  presented  in  this  report,  the  output  from  the  second  program  is 
used  directly  as  input  to  the  third.  However,  the  third  program  can  also 
accept  grids  generated  by  other  methods  provided  the  grid-control  functions 
are  evaluated  in  a  compatible  way. 

Calculations  have  been  performed  for  the  flow  over  the  stern  and  in  the 
wake  of  five  ship  hulls  for  which  experimental  data  are  available.  Extensive 
comparisons  have  been  made  with  the  data  to  understand  the  important  features 
of  such  flows  and  to  gage  the  success  with  which  these  can  be  predicted  by  the 
calculation  method.  The  results  are  presented  in  detail  so  that  they  can  be 
used  in  the  future  to  assess  other  methods  and  further  improvements. 

The  following  general  conclusions  can  be  drawn  from  this  study. 

(1)  The  present  numerical  method  has  many  attractive  features.  Among 
these  are  rapid  and  monotonic  convergence  to  steady-state  solutions  starting 
from  very  simple  initial  guesses,  ability  to  perform  accurate  calculations 
over  large  solution  domains  with  a  relatively  modest  grid,  and  the  potential 
capability  to  calculate  unsteady  flows. 

(2)  The  method  provides  an  excellent  description  of  the  overall  features 
of  the  flow,  including  the  pressure  and  friction  distributions  on  the  hull, 
and  the  mean  velocity  field  in  and  outside  the  viscous  layer  at  the  stern. 
Therefore,  it  can  be  used  in  practical  applications  to  calculate  viscous 
resistance,  determine  the  flow  ahead  of  propellers  and  appendages,  and  study 
the  problem  of  scaling  the  results  of  model  tests  to  full-scale  Reynolds  num¬ 
bers. 
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(3)  Calculations  for  a  typical  hull  take  about  5  to  io  minutes  on  a 
supercomputer.  It  is  also  possible  to  make  similar  calculations  on  a  modern 
minicomputer  in  a  matter  of  hours. 

(4)  The  basic  features  of  the  numerical  method  are  such  ’that,  with  minor 
modifications,  it  can  be  utilized  for  the  solution  of  many  other  types  of 
problems.  With  regard  to  the  numerical  content  of  the  method,  improvements 
can  be  made  in  several  areas,  the  principal  one  being  grid  generation  and 
control.. 

(5)  Improvements  in  the  description  of  the  flow  in  the  very  near  wake 
would  require  a  more  detailed  resolution  of  the  local  surface  geometry  and  a 
better  description  of  the  flow  in  the  near-wall  layers.  The  second  topic 
involves  further  developments  in  turbulence  models  for  near-waif  flows. 
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Fig.  9  Grid-dependence  tests  -  continued 
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Fig.  1  o  Influence  of  solution  domain  size  -  continued 
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Fig.  1 3  Girthwise  variation  of  pressure  and  friction  coefficients 
(i)  X  *  0.95;  (j)X«  0.9619 
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(a)  contours  of  axial  velocity  (U) 


Fig.  17  Velocity  field  in  transverse  sections:  measurements  of  Sarda  (1986) 

Top:  contours  of  axial  velocity  (U);  Bottom:  velocity  vectors  in  transverse  sections  (V,W) 
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Fig.  1 8  Profiles  of  velocity  components  (U,  V,  W),  pressure  (P),  and  turbulent  kinetic  energy  (k) 
compared  with  data  of  Sarda 
(a)  X  =  0.5 
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Fig.  19  SSPA  720  Cargo  Liner 

(a)  forebody  (b)  afterbody 
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Iscid  streamlines  and  measurement  stations,  Larsson  (1974) 
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Fig.  22  Girthwise  variation  of  pressure  and  friction  coefficients 
(g)  X  *  0.9;  (h)  X  -  0.9375 
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Velocity  field  in  transverse  sections ;  (g)  X  =  1 .0;  (h)  X  =  1 .1 
Top:  contours  of  axial  velocity  (U),  Bottom:  velocity  vectors  in  transverse  sections  (V,W) 
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Fig.  26  Contours  of  velocity  magnitude  from  experiments,  reconstructed  by  Broberg  and 
Larsson  <1984) 
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Fig.  27  Comparison  between  calculations  and  experiment 

(a)  Profiles  of  resultant  velocity  (Q)  and  crossflow  velocity  (C) :  X  =»  0.5 
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Fig.  27  Comparison  between  calculations  and  experiment 

(a)  Profiles  of  resultant  velocity  (Q)  and  crossflow  velocity  (C) :  X  *  0.75 
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Fig.  27  Comparison  between  calculations  and  experiment 

(a)  Profiles  of  resultant  velocity  (Q)  and  crossflow  velocity  (C) :  X  =  0.95 
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Fig.  27  Comparison  between  calculations  and  experiment 

(b)  Resultant  velocity  distribution  in  logarithmic  coordinates :  X  =  0.85 


Fig.  27  Comparison  between  calculations  and  experiment 

(b)  Resultant  velocity  distribution  in  logarithmic  coordinates  :  X  =  0.9 


158 


Fig.  27  Comparison  between  calculations  and  experiment 
(c)  Profiles  of  turbulent  kinetic  energy  (k) :  X  «  0.75 
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Fig.  29  Measurement  stations  in  the  experiments  of  Wieghardt  and  Kux  (distances  in  mm) 


Fig.  30  Typical  views  of  the  numerical  grid 
(b)  transverse  sections 
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Fig.  31  Pressure  and  friction  velocity  along  the  keel  and  waterline 
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Fig.  32  Girthwise  variation  of  pressure  and  friction  coefficients 
(a)  X=  0.6458;  (b)  X=0.7510; 
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Rg.  32  Girth'*/ise  variation  of  pressure  and  friction  coefficients 
(c)  X*  0.8745;  (d)  X=  0.9417 
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Pig.  35  Velocity  field  in  transverse  sections :  (b)  X  >  0.9103 

Top:  contours  of  axial  velocity  (U);  Bottom:  velocity  (V,W)  in  transverse  sections 
Left:  experiments  of  Wieghaidt  and  Kux;  Right:  present  calculations 
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Fig.  35  Velocity  field  in  transverse. sections  :  (c)  X  «  0.941 8 

Top:  contours  of  axial  velocity  (U);  Bottom:  velocity  (V,W)  in  transverse  sections 
Left:  experiments  of  Wieghardt  and  Kux;  Right:  present  calculations 
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Fig.  35  Velocity  field  in  transverse  sections :  (e)  X  -  0.9692 

Top:  contours  of  axial  velocity  (U);  Bottom:  velocity  (V.W)  in  transverse  sections 
Left,  experiments  of  Wieghardt  and  Kux;  Right:  present  calculations 
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velocity  and  pressure  fields:  profiles  of  U,V,W,  and  P 
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Fig.  36  Detailed  comparisons  of  the  velocity  and  pressure  fields:  profiles  of  U,V,W,  and  P 
(e)  X  -  0.9621 
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Fig.  36  Detailed  comparisons  of  the  velocity  and  pressure  fields:  profiles  of  U,V,W,  and  P 
(f)  X  *  0.9692 
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Fig.  36  Detailed  comparisons  of  the  velocity  and  pressure  fields:  profiles  of  U.V.W,  and  P 
(j)  X  =  1 .0440 
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Contours  of  longitudinal  component  of  vorticity;  experiments  (left),  calculations  (right) 
(b)  X  =  0.9493 


Contours  of  longitudinal  component  of  vorticity;  experiments  (left),  calculations  (right) 


Contours  of  magnitude  of  the  vorticity  vector;  experiments  (left),  calculations  (right) 


Fig.  38  Contours  of  magnitude  of  the  vorticity  vector;  experiments  (ieft),  calculations  (right) 
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Fig.  39  Offsets  of  SRI 07  Ore  Carrier 
(a)  forebody  (b)  afterbody 
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Fig.  40  Typical  views  of  the  numerical  grid 
(a)  symmetry  planes 
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Fig.  45  Velocity  field  in  transverse  sections  :  (a)  X  =  0.8 

Top:  contours  of  axial  velocity  (U);  Bottom:  velocity  (V,W)  in  transverse  sections 
Left:  experiments;  Right:  calculations 
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Fig.  45  Velocity  field  in  transverse  sections  :  (b)  X  =  0.9 

Top:  contours  of  axial  velocity  (U);  Bottom:  velocity  (V, W)  in  transverse  sections 
Left:  experiments;  Right:  calculations 
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Fig.  45  Velocity  field  in  transverse  sections :  (c)  X  *  0.95 

Top:  contours  of  axial  velocity  (U);  Bottom:  velocity  (V,W)  in  transverse  sections 
Left:  experiments;  Right:  calculations 
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Fig.  46  Detailed  comparisons  of  the  velocity  and  pressure  fields:  profiles  of  U,V,W,  and  P 
(a)  X  =  0.8 
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Fig.  46  Detailed  comparisons  of  the  velocity  and  pressure  fields:  profiles  of  U,V,W,  and  P 
(b)  X  =  0.85 
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Fig.  46  Detailed  comparisons  of  the  velocity  and  pressure  fields:  profiles  of  U,V,W,  and  P 
(e)X*  1.0 
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Contours  of  longitudinal  component  of  vorticity;  experiments  (left),  calculations  (right) 
(f)  X  =  1.0 


Fig.  48 


Contours  of  magnitude  of  the  vorticity  vector;  experiments  (left),  calculations  (right) 
(b)  X  =  0.9 


Fig.  48  Contours  of  magnitude  of  the  vorticity  vector;  experiments  (left),  calculations  (right) 
(c)  X  =  0.925 


Geometry  and  measurement  stations  in  the  experiments  of  Toda  et  ai.  (1988) 
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Fig.  51  Typical  views  of  the  numerical  grid 
(a)  symmetry  planes 
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Fig.  52  Pressure  and  friction  velocity  along  the  keel  and  waterline 


Fig.  53  Girthwise  variation  of  pressure  and  friction  coefficients 
(c)  X  =  0.95;  (d)  X  =  0.975 
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Fig.  56  Velocity  field  in  transverse  sections :  (a)  X  =  0.8 

Top:  contours  of  axial  velocity  (U);  Bottom:  velocity  (V,W)  in  transverse  sections 
Left:  experiments  of  Toda  et  al;  Right:  present  calculations 
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Fig.  56  Velocity  field  in  transverse  sections :  (c)  X  =  0.95 

Top:  contours  of  axial  velocity  (U);  Bottom:  velocity  (V,W)  in  transverse  sections 
Left:  experiments  of  Toda  et  al;  Right:  present  calculations 


2* 


T 


y 


y 


0.  fr - - 


.VANNS  v  vv  ^  ^  K 
— * ^\\VVv's\'s\NNS  V  Nw  vv 

^•.nNVnAaaVAW  WWvnnnv 
, .VVA\\\\\\\ \  \WS\Nnsw 
^-.^\\^\>\\\\\\\W  \  N  N  V  N  v 

/"s\  'W^'VNN' 


.  .  \ 
.  \  \ 
.  \ 


\\\\\\N\A\\* 

\M\\\\\\\\\ 
innwuu\ 

.  t  \  \  \  \  \  \  \  X  \  \  \ 

\  \  \  \  \  \  \  \  \  \  \  \  \ 

\  \  \  \  \  \  \  \  \  \  \  \  ' 

.  t  1  \  \  \  \  \  \  \  \  \  \  \ 

\  \  \  \  \  \  \  \  \  \  *  \  \  \ 

I  I  I  I  i  nu  m  n  « 

I!  M  u  n  u  »  m 

limmuM 

I'll  S80  <CAL. )  AT  X-.l 
III! 

tin  scale  «  uo-. 03 
4JlJLL_ - , — 


.  0« 


Y 


Fig.  56  Velocity  field  in  transverse  sections  :  (d)  X  *  0.975 

Top:  contours  of  axial  velocity  (U);  Bottom:  velocity  (V,W)  in  transverse  sections 
Left:  experiments  of  Toda  et  si;  Right:  present  calculations 
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Fig.  56  Velocity  field  in  transverse  sections  :  (g)  X  =  1 .00 

Top:  contours  of  axial  velocity  (U);  Bottom:  velocity  (V,W)  in  transverse  sections 
Left:  experiments  of  Toda  et  al;  Right:  present  calculations 


i 

! 


245 


* 


246 


Fig.  57  Detailed  comparisons  of  the  velocity  and  pressure  fields:  profiles  of  U,V,W,  and  P 
(b)  X  =  0.7 
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Fig.  57  Detailed  comparisons  of  the  velocity  and  pressure  fields:  profiles  of  U,V,W,  and  P 
(i)  X  =  1 .0 
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